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ON THE HIGHEST GRAVITY WAVES ON 
DEEP WATER 


By HAROLD JEFFREYS (St. John’s College, Cambridge) 
[Received 2 November 1950] 


SUMMARY 

The method of Michell and Havelock in discussing the form of the highest possible 
progressive wave on deep water depends on the rearrangement of a double series and 
the determination of a set of coefficients by successive approximation. At one point 
in the method a numerical check gave a discordance. It is found here that the 
discordance is probably due to the retention of a small number of terms in a slowly 
convergent series; it disappears when the limit is estimated from the known 
approximations by L. F. Richardson’s method of deferred approach to the limit. 


By a well-known result of Stokes, if a sharp crest forms on a gravity wave, 
the angle at the crest should be 120°. The presumption is that with 
increasing ratio of amplitude to wave-length the crests become sharper and 
the troughs flatter, and that there is an extreme possible height such that 
the crest becomes angular. With additional supply of energy the wave will 
be unable to store the energy, and this is associated with the formation of 
‘white horses’. Stokes and Rayleigh have shown, by expansion in powers 
of a function of the amplitude, that the waves do behave in the way 
indicated for sufficiently small amplitude; a rigorous proof has been given 
by Levi-Civita, but it is still unknown up to what amplitude the expansions 
are valid. 

J. H. Michell (1) attacked directly the problem of gravity waves when 
the crests are already sharp, and derived an estimate of the maximum 
ratio of height to amplitude. His method has been extended by T. H. 
Havelock (2), who has also made some numerical corrections. The maxi- 
mum ratio is found to be about 1/7. The result has, however, been queried 
for two reasons. From observation the critical ratio appears to be substan- 
tially less, values near 1/12 being quoted. From theory Michell and Have- 


lock give two methods for estimating the departure of the slope from 


uniformity near a crest and on one side of it, and the results do not agree. 
This failure of a check might suggest an error in the work at some point. 
The object of this note is to point out that the disagreement is due to the 
slowness of convergence of one of the series, and that when this is allowed 
for by L. F. Richardson’s method of ‘deferred approach to the limit’ the 
agreement is as good as could be expected. 

[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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Michell’s method is one of conformal representation, in which the 
variable z is expressed in terms of the complex variable w = ¢-++-is; the 
free surface is 4 = 0 and the liquid occupies the region ys > 0. The prob- 
lem is re-scaled so that the wave-length for 4 is 7. He finds 


dw eo ve vee oe 

; (—isinw)te"(] 6,e"* ! b,c“ a (1) 
dz 

which has the period z in the region 0. The b, are real. He gives the 

] ’ 


condition at the free surface as 
49 —, (2) 


where q is the speed of the liquid. In this equation each side is expressed 
in the form sin‘¢ times a series of cosines of odd multiples of 4, which 
therefore appears to preserve the correct periodicity. The coefficients b, 
are determined by equating coefficients and solving numerically. The 
reduction of the left side of (2) to the required form is straightforward but 


complicated, since the coefficients are of degree 4 in the b,. But 


c d 


-sin'd{sin(4¢- 17) > b,.sin(2rd : ld Lor), (3) 
cy 


Ly 
which, if sin'd is taken to be real, has period 37 instead of 7. In fact the 
taking of the imaginary part of (1) over % = 0 requires attention to the 
behaviour of sin'w near its singularities, and (3) is correct only for 
0] b< 7. 

Michell’s next step is to expand each sine in the series in (3) in a cosine 
series for 0 < 4 < 7; by rearrangement he finds the series actually used. 
As values of ¢ outside the interval (0,7) are not used the ambiguity of (3) 
outside this interval does not affect the expansion and his method should 
be correct subject to convergence conditions. 

However, the coefficients of cos(2n--1)¢é derived from each term of (3) 
decrease with n like n~*, and it may be expected that the error arising in 
the solution through taking n terms at ¢ = 0 would be of order n~'; this 
error is the result of the rearrangement of a double series and has nothing 
to do with the accuracy of the substitution of (1) into (2). In fact Michell’s 
method for finding the height of the wave (crests above troughs) depends 
on the series 1—b,+b,—b,+..., and we should expect the error due to 
stopping at the nth term of this series to be of order n~? (since the b, are 
all positive). Of his two methods for studying conditions near the crest, 
one should give similar convergence and gives ¥ b, = 0-0758. The other, 
using the numerical values of the 6, leads to ¥ b,, = 0-0511 from three 
terms of the series. 
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GRAVITY WAVES ON DEEP WATER 387 
Havelock extends the analysis to waves of less than the critical height and 
ympares the results with those of Stokes’s development. The convergence 
cessive approximations used in calculating the maximum height 


ppears to be satisfactory and Michell’s ratio 0-142 for height /wave-length 


is confit d. But he finds that ¥ 6, increases with successive approxima- 
tion. Keeping only 6,, he finds b, = 0-0311. When three terms are kept,t+ 
0-0540: and for four terms. b, bh, bs b, 0-0584. Now 
suppose that the sumS,, of the series b, +-b,+-b,+-..., when n terms are kept, 
' S A Bn-1+Cn-?. 
Fitting this form to the three values found gives 
N 0-0737—0-0674/n+0-0248 n2. 


‘he limit 0-0737 obtained by this method of ‘deferred approach’ agrees 
well enough with Michell’s value 0-0758, considering the roughness of the 
extrapol 


\ccordingly, there appears to be no reason to doubt the theoretical 


ilidity of the Michell—-Havelock method, and, as it coalesces with the 


Stokes method for the small amplitudes where the latter is certainly valid, 
ere is a strong presumption that it is correct over the whole range of 
possibl litudes. The explanation of the disagreement between theory 


ind observation is likely to be found elsewhere. 


Since the slowness of the approximation consists in changes of the indi 
lual b. at successive approximations, the 5, themselves in each approxi 


decreasing rather rapidly with 7, it appears that a more direct 


hod would have been better. It would be possible to assume a finite 
umber of tern in (1), and to choose the coefficients to fit (2) at a finite 
number of points; g then emerges as an eigenvalue. The difficulty of the 


solution would come from the fact that q* is of degree 4 in the b,, but this 
is already present in Michell’s method and could probably be taken into 
vecount at st as easily in a numerical solution. It is possible that the 
successive approximations to the maximum height itself are not yet 


suincient! 


alues did not completely satisfy his equations, 
{ I Havelock’s values. Havelock’s B, are Michell’s by. 


REFERENCES 
i. J. H. Mi Phil. M 36 (1893), 430-7. 
A, 95 (1919), 38-51. 














ON VIRTUAL MASS AND TRANSIENT MOTION IN 
SUBSONIC COMPRESSIBLE FLOW 


By JOHN W. MILES 


(Department of Engineering, University of California, Los Angeles) 
[Received 30 January 1951] 


SUMMARY 

The transient motion of a body in a compressible fluid is discussed on the basis of 
the acoustic approximation, and it is shown that the concept of virtual mass is 
applicable only to steady flow. It is found that if a velocity is imparted suddenly 
to a body half the work is lost in radiation. The velocity and acceleration resulting 
from a constant force suddenly applied to a circular cylinder are calculated with the 
aid of the Laplace transformation, the results being given in the form of curves and 
also as expansions in ascending and descending powers of t. The related case where 
a force is applied suddenly at ¢ 0 and removed at ¢ t, is treated, and the radia- 
tion loss is calculated as a function of ¢,. Also treated is the case of a suddenly 
imparted velocity, which is an extension of the classical discussion of impulsive 
motion. It is concluded that the radiation loss due to acceleration is not likely to 
be important in practice, with the possible exceptions of airships and submarines. 


Notation 
F(s) Laplace transform of f. 
I impulse. 
I, Bessel function of imaginary argument and order one. 
K, Macdonald function of order one. 
M Mach number. 
P force per unit length applied to cylinder. 
constant force. 
V(s) Laplace transform of v(7). 
W work (per unit length) done in accelerating body. 
Y(s), Z(s) see equation (2.11). 
a radius of circular cylinder. 
c sonic velocity. 
f(z) dimensionless force; see equation (2.7). 
m virtual mass (per unit length). 
p perturbation pressure. 
Po Static pressure. 
r polar radius. 
s Laplace transform variable. 
t time. 
t, time at which force ceases to act on cylinder. 


{Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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wu velocity of cylinder (normal to axis). 
constant velocit 
u final velocity. 
acceleration. 
dimensionless velocity. 
[ gamma function 
Y'(s) Laplace transform of u(r). 
ratio of specific heat at constant pressure to specific heat at con- 


stant volume. 


starting efficiency (ratio of final kinetic energy of flow to work 

53 done). 
lenly 6 polar angle. 
Iti . € dimensionless, polar radius. 
p mass density of fluid. 
ne o positive, real constant. 

7 dimensionless time. 
er d velocity potential. 


y% dimensionless velocity potential; see equation (2.4). 





es. 1. Introduction 

Tue problem of impulsive motion of a body in an ideal, incompressible 
fluid has been discussed by Lamb, Larmor, Kelvin, and others (1). As 
| might be expected from the concept of virtual mass, an impulsive force is 
required to generate impulsive motion, in so far as the fluid is assumed to 
be incompressible. However, if the compressibility of the fluid is taken 
into account the force required to generate impulsive motion may be ex- 
pected to remain bounded and to vanish only asymptotically with time. 
Moreover, the concept of virtual mass, while still valid for the description 
of the energy and momentum of the fluid flow in the steady state, is of no 
iid in computing the transient response. In particular, we shall see that 
the work required to generate the impulsive motion is not entirely recover- 
tble, so that the process is irreversible, and some part of the starting force 

must be regarded as ‘drag’.’ 
To illustrate these phenomena, we shall consider here the transient 
motion of a circular cylinder at low, subsonic speeds in an ideal, com- 
pressible fluid and shall treat in detail the particular cases of a suddenly 


applied force and a suddenly imparted velocity. (The solutions to these 


i 
We remark that this paper was motivated by a brief discussion of just this question 
Prof n Karman, who raised the question of the validity of the virtual mass concept 

4 Supersor flow. Our thank are due to Prof. H. K. Forster, with whom we dis- 

ed the sami iestion. Prof rster pointed out the similarity between the present 
celeration of an electron, with which is associated radiation. We are 


nee 7. 
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problems also may be of interest in acoustics.) Actually the force acting 
to produce impulsive motion is simply the Bromwich integral of the force 
acting on a harmonically vibrating, circular wire, the latter problem having 
been treated by Rayleigh (2). However, the evaluation of this integral is 
rather involved; indeed, the difficulties presented by this evaluation have 
rather discouraged us from attempting a solution to the associated problem 
of the scattering of a sharp-edged plane wave by a circular cylinder. In 
this connexion we remark that Fox (3) has treated the analogous problem 
of scattering by (or impulsive motion of) an infinite plane ribbon, using an 
iteration procedure similar to that proposed by Schwarzschild (4) for the 
harmonic case. Unfortunately this procedure does not lead to a solution 
of the latter (harmonic) problem in closed form, a better approach being 


offered through separation of the wave equation in elliptic coordinates (5), 


and it does not appear that our present knowledge of the properties of 


Mathieu functions would allow any direct evaluation of the appropriate 
Bromwich integral. 

With reference to more general transient problems, we remark that their 
solutions may be synthesized from the ‘step function’ solutions with the 
aid of the Duhamel superposition theorem. 

In the case of aircraft, the accelerations actually experienced in subsonic 
flow are likely to be so small as to render the considerations discussed 
herein of academic interest only. While acceleration effects might be 
relatively more important for airships or submarines, it appears that the 
most likely practical applications of the results to be obtained would be 
in the field of acoustics, where the linearization of the equations is entirely 
appropriate. 

Before attempting detailed formulations we examine some of the simpler 
aspects of the impulsive motion problem which may be established with a 
minimum of analysis. First, after equilibrium is attained the flow will be 
potential, and there will be no resultant force on the body. Secondly, if m 
is the virtual mass of the body in this potential flow, the final momentum 
of the associated flow will be muy. where uw, is the impulsively imparted 
velocity; hence the total impulse delivered to the body by the transient 


force, P(t), required to impart the motion is given by 


l ( P(t) dt = muy. (1.1) 


Thirdly, since this force acts at a constant velocity, the total work done 


is given by 


W = [ P(t)uy dt = mu2. (1.2) 
! 0 ( 


0 
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Now the energy of the associated flow after the steady state has been 
; whieved is }muj. so that an equal amount of energy evidently is dissipated 
during the period of transient motion. The agency of dissipation is, of 


eours radiation Mor 


kinetic energy smus as 


somewhat 


ent problem of 


‘over, if the wire were to be stopped suddenly, the 
We 


ilogous situation exists in the rather more 


d with the flow also would be radiated. 


r system of one degree of freedom. For ex- 


nea 
aati 


mple, if an ideal capacit is charged by a suddenly applied d.c. source, 
= this source delivers exactly twice the final energy stored in the capacitor, 
the even though there be no hmic) resistance in the circuit, the energy 
tio difference being accounted DY radiation. 
= We rem that the results (1.1) and (1.2) are valid only in so far as up 
5 is small compared to th mic velocity, c, viz. 
ol M Uy /C 1. (1.3) 
put are itherwise indepet ient ot ¢. 
heit 2. Transform solution of problem 
the The behaviour of small disturbances in a frictionless fluid is governed 
by the wave equation 
VY" dy hd A(r. 0.1). (2.1) 
where the velocity potential, c the sonic velocity, and (r,@) polar 
th coordinates, as shown in Fig. 1. If py and p, are the pressure and mass 
density of the fluid at rest and y the ratio of specific heats, ¢ is given by 
c* = YPo! Po: (2.2) 
nd the 7 irhation pre ( \ 
al : 
e 
p Po M- (2.3) 
7 n ordet render our problem dimensionless, we choose the radius of 
he cireul vlinder. a S ‘haracteristic length and ¢ as a characteristic 
nn velocity. with (a/c) result a characteristic time. Moreover. since we 
badd shall deal only with motion (of the body) in a fixed direction (0 0). we 
! may suppress the angular dependence (cos @). Thus, we introduce a modi- 
ed, dimensionless potenti such that 
P(r, 0,1 a,ct,a)cos 0, ub h(E, 7). (2.4) 
LI Substituting in (2.1) vields 
us__) Eu us 0). (2.5) 


condit 


} } | ‘ 
HOt aary 


prescribed velocity. u(t 


latter 


resulting 


irom a 





] 7] 
LUI 


it the circular cylinder is dictated either by a 


or a prescribed force (per unit length), P(t), the 


Introducing the 


peripheral integration of (2.3). 
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dimensionless quantities 7(7) and f(z), we may write these boundary con- The | 
ditions in the form - 
bel(E,7)\e-4 u(r) u(t)/e, (2.6) first 
. while 
(1,7) = f(r) = P(t)/zaypyp. (2.7) ; 
force 


In addition to these conditions, we require the solution at infinity to behave 
as an outgoing, cylindrical wave (the Sommerfeld radiation condition). 
Our problem is then reduced to relating v(7) and f(z). 





Fic. 1. Cross-section of circular cy linder, showing polar coordinates. 


A formal solution to the foregoing problem can be obtained directly, 
using the Laplace transformation. Thus, writing 
¥'(é, 8) Pib(E, 7)! | e*™h(E, 7) dr (2.8) 
. 3.7 
W(é, 8) [V(s)/sK4(s)]K,(s€) | F(s)/sK,(s)|Ky(s€), (2.9) L 


where K, is Macdonald’s solution to Bessel’s equation of order one, and 


we obtain 





V(s) and F(s) are the Laplace transforms of v(r) and f(z) respectively. 
Then, eliminating K,(sé), we have 








V(s) = sY(s)F(s), (2.10) whe 

sY(s) | K}(s) K,(s) ] sZ(s), (2.19%) real 

where sY(s) and sZ(s) are the ‘admittance’ and ‘impedance’ functions this 
respectively. The formal solution is completed by applying the Faltung the 
theorem to (2.10) to obtain in tl 
: the 

v(r) = f(O+)y(r)+ | y(r—A) f(A) dA, (2.12) aia 

. : as t 

f(r) = v(O+)ze(r)+ | 2(7—A)v’(A) dd. (2.13) give 


and 
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on The function y(7) Y-\‘)'(s)} is essentially the ‘indicial admittance’ of the 
T first problem, i.e. the (velocity) response to a suddenly applied unit force, 
i while 2(r) has a similar interpretation for the second problem, being the 


force response to a suddenly imposed unit velocity. 
; 


- Im(S) 


oo 
a 








ae a RAS) 








thy 
Fic. 2. The complex s-plane, cut along the negative real 
xis, showing the path of integration for the various in- 
. verse transforms. 
y. sS 
3. The evaluation of y(r7) 
2.9 Taking the inverse transform of Y(s), as given by (2.11), we have 
and ix 
ely y(7) | A4(s) sK,(s) |e’? ds, (3.1) 
10 ‘ ‘ : 
where o is positive, real constant, and the s-plane is cut along the negative 
ll real axis, in recognition of the logarithmic branch point in A,(s). To reduce 
ions this integral, we consider the contour integration of the integrand around 
ung the path ABCDEFGH, shown in Fig. 2. Now, since K,(s) has no zeros 
n the finite, complex plane (6), the integrand is regular everywhere within 
19 the contour, and the complete, contour integral vanishes. Moreover, it is 
12) , 


| readily shown that the contributions of the paths BCD and GHA vanish 
is the contour is allowed to expand indefinitely to the left. Hence y(t) is 
13) given by the negative sum of the contributions of the paths DE, FG, 


ind EF, 
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The integral around FE yields a term 7. The integrations over the paths | and 
DE and FG can be simplified with the aid of the result (6) 
Kk, (ue=7") Ky (wail (a), (3.2) 
where J,(u) is Bessel’s function of order one and imaginary argument, and wher 
by using also the Wronskian relation between A, and J, (6). The final | 
result is 
Yy(T) 7+ u . KF(u) a TF(u)| lp—uT dy, (Ss. 
0 
lO 
8 
7 
— 
Es 
“T 
O ‘ 
Fic. 3. y(7)—7 as given by equation (3.3). =e 
expa 
The integral in (3.3) apparently cannot be evaluated in terms of simple | totic 
known functions. However, precisely the same integral has been evaluated 
by the British Admiralty (7) in connexion with Lighthill’s analysis (8) of y 
supersonic flow past a ‘quasi-cylinder’, being there designated as l’(2). 
(The actual evaluation was carried out by numerical integration of an 
associated integral equation.) wher 
The function y(7)—7, as given by (3.3), is plotted in Fig. 3 for values | lave 
of 7 between 0 and 10. Also plotted, in Fig. 4, is y'(7), given by logar 
. 
y (tr) l | wl AG (w) +a? TF (uw) |te-" du. (3.4 
0 
The physical solution represented by Figs. 3 and 4 corresponds to the with 


sudden application of a force, P,, at = 0 with a resulting velocity 


u(t) (Pye TAY Po )y(ct/ a) (vu 
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nd an acceleratiol , ean 
ar — tA m)y (ct/a), (3.6) 
7a" p, (3.4) 
cher is the apparent mass in the final potential flow. We note that 
ce y' (0 the initial apparent mass is twice the final value. 
/ 
8} / 
| / 
| / 
Gi 
ee | 
= § 
> 4| 
vi 7 
| 
; > : 
O 4 6 8 lO 
‘IG. 4 is given by equation (3.4). 
‘or some purposes it may be convenient to have series representations 
(7). Appropriate expansions for small and large 7 can be obtained by 
xpanding ) for large and small s, respectively. Introducing the asymp- 
mp totic and ascending expansions of A,(s) in (2.11), we obtain 
) s1+4s-14 352+ O(s-3)], (3.8) 
8 , 
; Y(s C'+-In(s/2 s*| (2C2—2C+1) 
it 2(2C’— 1)In(s/2)+-2 In?(s/2)|4-O(s* Ins), (3.9) 
e C is Euler’s constant (0-5772...). The inverse transforms of the 
erse powers of s are elementary (9) while the inverse transforms of the 
ithmic terms can be obtained by differentiating the result (9). 





sin(7a) ; 
) ot = (@ Ll), (3.10) 
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4. The evaluation of z(r7) 5. I 
From (2.11) we have I 
i - the 

os 27 [ Ky (s)/sK4(s) Je? ds. (4.1) |} furt 

o—i« tion 


Again, we consider a contour integration around the closed path of Fig. 2 
the integrand of (4.1) being regular in the cut plane of that figure. In this 
case the integrand is regular at the origin, so that the contribution of the 
path EF is null. On the other hand, A}4(s) has two zeros in the left half- 


plane, which may be determined from the known zeros (10) of H(z) to | 
be at the conjugate points 
ifs, ) = 0; &=—% 0-6435-+7 00-5012. (4.2 
The corresponding residues, simplified by the use of the differential equa- 
tion and recurrence formula for A,(s), are found to be 
Res | — Ky(s)/sK4(s) |,—<, | Ay (8,)/s, K4(s,)} —s,(1+s?)-!. (4.3) 
The contributions of the paths DE and FG can be evaluated as in the | 
previous case, and the end result is 
2(7) exp( 0-64357)| *2120 cos(0-50127)-4 0-1898 sin(0-50127) | : 
— | u*[ KP(u) +a lP(u)|-te-"* du. (4.4 
0 
This function, too, was calculated in (7), being designated therein as f(a 
and found to be a solution of the integral equation at 
: are 
| k(r A)z(A) dxr [r(7 +2) ]?, (4.5a) 
0 
. ’ ] 
where k(r) (r+1)*[7(7+2)]-?. (4.5b 
fro 
The result is plotted in Fig. 5. how 
Interpreted physically, Fig. 5 gives the force-time history required to 
impart suddenly a velocity uw, to a circular cylinder in accordance with 
the relation Ex 
P(t) = raypo(Up/c)z(ct/a). (4.6) 
It is of particular note that this force reverses at about 7 3-13 and we 


thereafter remains negative, although approaching zero asymptotically. 
To obtain expansions for small and large 7, we proceed as in the previous 
case, finding Pg Ni ‘ 
1 2(r) = 1—}r— ger? + fy? + O(r4), (4. 


(4.5) 
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5, Rectangular force pulse 
[t was shown in the introduction that the impulsive starting process of 
the preceding section is only 50 per cent. efficient. In order to obtain some 
+. further correlation between this starting efficiency and the rate of accelera- 


tion. we consider the case where a force P, is applied to a circular cylinder 
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Fic. 5. 2(7), as given by equation (4.4). 
us 
tt = 0 and removed at ¢,. The corresponding velocity and acceleration 
ire given by . 
“s : u(t) (P, c/mayp,)| y(ct/a)—y(e(t—t,)/a)], (5.1) 
1.5 acct , Se 
u(t) (P,/m)|y'(ct/a)—y'(e(t—t,)/a)]. (5.2) 
is From direct consideration of the total impulse delivered, or otherwise 
from an examination of the Laplace transform of (5.1) in the neighbour- 
hood of s 0, the final velocity is given by 
Da ‘ — 
— u lim w(t) = Pyt,/m. (5.3) 
wit itis 
Examining the first two terms (for small s) in the Laplace transform 
1 6 
fSy(r)—y(t—7)} = [1 —exp(—sz7,)|¥(s) (5.4) 
we find that this final velocity is approached in accordance with 
lly 
U(r) P " —— 
evious y(t l—7+-° O| (r T}) MIn(7 7,)}. (9.9) 
} The final kinetic energy associated with the potential flow is given by 


4.8 Smus (yt,)" 2m. (5.6) 
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On the other hand, the work done by P, is given by 


ty 

W | P, u(t) dt n—1(mu2,/2), (5.7) 
0 
»n = 473/ | y(t) dr, (5.8) 


where 7 is the starting efficiency and is plotted in Fig. 6. For small values 
of 7,, we may integrate (3.11) directly to obtain 





7) AT, Lert O(74). (5.9) 
— — 
st —_—— 
| | 
6 Sa ee J 
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| = 
O . * 6 8 10 
Ty 
Fic. 6. (71), as given by equation (5.8). 
For large values of 7,, we have 
TL {y-t} = 28-1¥(s). (5.20 


> 


Substituting Y(s) from (3.9), we may interpret the various terms as in that 
section, excepting s-!Ins, which, however, is a tabulated transform (9 
Thus, we obtain 
” 1 — 277? In(27,)+77 4[ 4 1n?(27,)+ 2 In(27,)—1]+ O(77* In?7,). (5.11 
[t is evident from (5.9), and also from Fig. 6, that the starting efficiency 
of a short pulse, entailing as it does high acceleration (the initial accelera- 
tion being w, 2t,). is quite poor, while for long pulses (large 7,) the efficiency 
approaches unity. However, we note that, even for large 7,, there is always 
a finite amount of energy dissipated, since 
W- imu? 


(P5/7ypo)|n(2ct,/a)+ O(77 4 In?7,). (5.12 
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The limiting case of the rectangular pulse under discussion leads to the 


impulsive force, defined by 


lim | P(t) dt = Pt, = mu,. (5.13) 


The resulting velocity, as a function of time, is obtained by taking the 


\pproprl te limit in (5.1 vnd leads to 
u(t) i 7 
y (ct/a), (5.14) 


u 


which has been plotted in Fig. 4. For small time this result yields 


(ft) 


1+ 874 O(7?), (5.15) 


while for large 7 it is identical with the result for a pulse of finite duration. 
We emphasize that this type of impulsive loading is not that contemplated 
hv Lamb ef al.. which is rather the limiting case of section 4 as c is allowed 
to become infinite. In the latter case the initial velocity is, of course, 
identical with the final velocity, and the efficiency is 50 per cent. 

In order to form some estimate as to the importance of the effects dis- 
ussed herein in typical physical problems, we consider approximate orders 
f magnitude for 7,. The order of magnitude of ¢, is usually the same as 

of so that we may write 


O}| (c?/ag)(u/e)(g w)|. (5.16) 


For an aeroplane accelerated along the line of flight (to which, strictly 
speaking, our results are not applicable because of the linearization) the 


l 


characteristic length should be taken as the wing thickness or the product 
of the chord and the angle of attack, either of which might be about 1 tc 
2 feet for typical configurations, and c?/ag is then in the neighbourhood of 
>to 4x 104. Since (u/c) and (g/v%) might be of the order of 0-1 and 10, 


respectively, in horizontal flight, or 1 and 1 in a dive, it is evident that 7, 


will be so large as to give | for all practical purposes. On the other 


hand. in the ease of a large balloon undergoing its initial acceleration, we 
might obtain values of 7, as low as 4, which would correspond to a starting 
efficiency in the neighbourhood of 75 per cent. Turning to an acoustical 
example, ¢, might be of the order of 10-3 sec. for a violin string, with a 


bout 10-4 feet, so that +, would be about 104, while for a large drum 7, 


6. Conclusions 
We emphasize that the practical significance of the foregoing calcula- 


tions is limited by the initial linearization, which neglects second-order 
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terms in velocity. Nevertheless it may be expected that they yield a 
picture that is qualitatively correct throughout the subsonic régime, and Al 
we conclude that, while the acceleration of a body at low speeds in a. per- 
fect fluid is to some extent irreversible, the fraction of the expended energy 
rendered unavailable is not likely to be of importance in most practical \ 





problems, with airships and submarines as possible exceptions. 


(Added in proof. Prof. Lighthill has informed us that one of his 
students has carried out a similar analysis for a sphere.) 
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and |} AN EXPANSION FORMULA FOR THE DRAG ON A 
i CIRCULAR CYLINDER MOVING THROUGH A 
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fi VISCOUS FLUID AT SMALL REYNOLDS NUMBERS 
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xj By S. TOMOTIKA and T. AOI 
) 
(University of Kyoto, Japan) 
Received 30 January 1951] 
SUMMARY 
In the present paper an expansion formula is derived for the drag experienced by 
ular cylinder moving through a viscous fluid at small Reynolds numbers, by 
ting from the general expression for the drag on the cylinder obtained by us 
isly. Our expansion formula is expressed in a power series of the Reynolds 
ACS number of the undisturbed flow and includes Lamb’s well-known formula as its first 
approximation. The range of values of the Reynolds number for which our expansion 
, Brit formula is applicable is estimated. We find that our expansion formula correct to 
fourth power of the Reynolds number can be used with sufficient accuracy, 
( provided that the Reynolds number is numerically less than 4. 


1. Introduction 

Tue steady flow of a viscous fluid past an obstacle can be discussed satis- 
factorily on the basis of Oseen’s linearized equations of motion, provided 
that the Reynolds number is fairly small. Recently (1,2) we have obtained 
an exact solution of Oseen’s equations for the case of the steady flow of an 
incompressible viscous fluid past a circular cylinder, and the flow patterns 
round the cylinder have been computed. We have obtained a general 
expression for the drag experienced by the cylinder, with special reference 
to the pressure drag and the frictional drag. The values of the drag coeffi- 
cient for several small Reynolds numbers have been computed numerically 
and it has been found that the calculated values are in good agreement 
with observation. In numerical computation of the values of the drag 
coeflicient, however, a system of simultaneous linear algebraic equations 
has to be solved numerically for each given value of the Reynolds number 
so that the work is rather cumbersome. It is therefore desirable to derive 
in expansion in powers of the Reynolds number, which is convenient for 
numerical computation of the values of the drag coefficient. The object 
of the present paper is to derive such an expansion. As is well known, the 
corresponding formula of expansion for the drag on a sphere was obtained 
by Goldstein as early as 1929 (3). 

[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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We here obtain a formula of expansion correct to the fourth power of the 
Reynolds number for the drag coefficient of a circular cylinder. As would 
be expected, Lamb’s well-known approximate formula is included in our 
formula as the first approximation. 

Recently a similar investigation on the drag on a circular cylinder has 
been made by Sidrak (4), who has also derived a formula of expansion for 
the drag coefficient. But on examination we find that his result is 


erroneous. 


2. An expansion formula for the drag 

For convenience we first refer to the general expression for the drag on 
a circular cylinder, obtained in our previous paper (1). We assume that 
a circular cylinder of diameter d moves with constant velocity UU in an 
incompressible viscous fluid with the kinematic coefficient of viscosity y, 
and we define, as usual, the Reynolds number R of the fluid flow as 
R = Ud/v. The drag experienced by the cylinder is then given by 


> y % 2 
D 2rpyl > B..., (1) 
m-0 
where , is the coefficient of viscosity of the fluid. The constants B,,, are 


determined by solving the following simultaneous linear algebraic equations 


¥ Badna(B) =| fy — 2, 2 


m=—0 (n ao Sn 2 


The coefficients A,,,, are functions of the Reynolds number F and are 
expressed as 
r Dn —nKm—1t nen Kmaittn—nis Amt ,K 


mn mn m man m m-n m m+n m? 


the argument of the modified Bessel functions J,, and K,, being } R. 


As usual we introduce the drag coefficient C), defined as C;, = D/(pl*d), 
p being the density of the fluid. Then, the general formula for C,, becomes 


5 a 
’ «71 ) . 
Cp : R S B... (9 
UC 
m=0 


In order to calculate the drag coefficient for any given value of the 
Reynolds number, we have first to determine the B,, by solving the 
above system of simultaneous equations (2). Theoretically the solution 0! 
(2) can be obtained by means of infinite determinants; in practice, to find 
the numerical solution for any given value of the Reynolds number wé 
solve a finite number of equations. Thus, for the first approximation, we 
put B,, By, Bs,... equal to zero, and solve the first equation for By. For 
the second approximation we put B,, Bg,... equal to zero, and solve th 


first two equations for B, and B,, and so on. The first, second, and third 
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approximate values of the drag coefficient plotted against the Reynolds 
number have been given for each approximation.t 

For our present purpose, that of deriving an expansion formula for the 


drag starting from the general expression (1), or (3), we first express the 


\.. in power series of R by using the expansions of the modified Bessel 
functions I and kK, - Denoting for brevity 
S ‘—v—log 1 R = 2-0022—log R, 

where 7 0-57721... is Euler’s constant, we obtain the following expan- 
sions P 1) D2 7 . oe 4 

Nor = 28+ (gs +dsS) R2+(alss-t adhe S) RY+-...; 

l 1YV » 1 1 y 23 
) (—744+48)R+(25+-3,S) R34 
\o ay) R? 
1 22 
Ait i+ S)+ 92 
l > 
A, ) 4 : UR 
2 R’? 

\ l 

si 2 

As 

l 
hoe 8 
= R 
| 
i. «= 
7 hk? 


\ssuming R to be small, it is easily seen that the orders of magnitude of Bo, 
B,, B,,... ave 1, R?, R*,..., respectively, and that the method of approxima- 
tion mentioned above can be applied appropriately to our present calcula- 
tion. Making use of the expansions for the A,,,, just given, expansion 
formulae for the drag coefficient corresponding to the first, second, and 
third approximations can be obtained. 
The first approximation gives 
477 
Cp RS’ (4) 


and this is the well-known formula obtained first by Lamb. 


Nl: Table I and Fig. 5 in ir paper (1). 
ig l 1 ilat the series definitions for J,, and K,, are those adopted in 
G. N. Wat I The j of Bess Functions (Cambridge, 192: 
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Next, the second approximation is expressed as 


; 4 De ia ge ee ‘ 
Cp RS ae 3S . i) 3a}: (5d) 





Further, the third approximation yields an expansion formula of the 


form? 


4a f P ? ] Ge  . 
Co _ (S?—4S8-+ +) R (S#—18384+ 2S— 3) . . (6) 
RS S 5 ee : = midis Yn 





The fourth and higher approximate formulae may be obtained in like 
manner, but it will be seen in the subsequent section that the above third 
approximate formula (6) is sufficient for computing values of the drag 
coefficient when the Reynolds number is fairly small. 


3. Numerical computations 

It is worth while estimating the range of values of the Reynolds number 
for which our expansion formula (6) is applicable. This can be done by 
examining the convergence of (6) as well as by comparing the values of the 
drag coefficient calculated from (6) with those calculated from the general 
expression (3) by solving the system of simultaneous equations (2) numeri- 
cally for B,, B,, and B,. In order to examine the convergence of formula 
(6), the values of the three terms, corresponding to the orders of magnitude 
of 1, R?, and R* respectively, in the square bracket in (6) have been 
calculated numerically for various values of the Reynolds number less 
than 6. No numerical calculation has been made for values of R greater 
than 7, however, since the quantity S vanishes at R — 7-4 approximately. 
The results are given in Table 1, from which it will be seen that our 
formula (6) begins to diverge at R 4, below which, however, it con- 
verges rather rapidly. 

TABLE | 


R rst term 2nd term 3rd term 


“4 I I 000 
*( I 2 0o*00 
8 I 0°04 0°00 

I I 05 0°00 
I oO'l3 o"02 

I O'2I1 0°04 

| I O°31 0°02 
5 I O54 0°22 
6 I 1°34 pr 2°24 


The values of the drag coefficient calculated respectively by the first. 
second, and third of the approximate formulae ((4), (5), (6)) are given in 
+ We have recently obtained the corresponding formula for the drag on an elliptic cylinder 
whose major axis is either parallel or perpendicular to the undisturbed uniform flow, and we 
have confirmed that the above formula (6) agrees, as would be expected, with a limiting 
form of the corresponding expansion formula for the drag on an elliptie cylinder. 
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Table 2 and are shown graphically in Fig. 1. 
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The last column of Table 2 


























~ gives for comparison the values of the drag coefficient computed numeri- 
cally by the general formula (3). It will be seen that the values of C), 
the calculated by the third approximate formula (6) are in good agreement 
E = 
6 Cp 
\ | | 
ro 10 \ + ———E—EE + | Bcnsmmscmnemenel 
like } 
urd \ | 
1 & \ _ ——+ —_—__—— 4 
\ | | 
6 
l 
1D¢ 
i 
4 
TN 
he 
el 
1e1 a : ae, a ae | 
nu | | | 
hail 0 | eee Eee | il 4 
0 | 2 3 4 
beet R 
leas Fic. 1. Curves of the drag coefficient Cp, plotted against the Reynolds number R. 
; Curves I, II, and III represent respectively the calculated values of Cp by the first, 
l l 1 ~ > . . ° 
md, and third approximate formulae ((4), (5), (6)), while the remaining curve 
U shows the values of Cp calculated by the general formula (3). 
‘i 
with those calculated numerically by the general formula (3) in the range 
of values of the Reynolds number in which our expansion formula (6) con- 
verges rapidly. Thus, we conclude that our expansion formula (6), correct 
TABLE 2 
Values of ¢ D 
1) (2) 
10°63 I 3 
“I “13 13 
g 65 6:78 
$07 $04 
first 547 39 
en 1 ai abo 
ore as far as terms of the order of R*, gives values of the drag coefficient with 


sufficient accuracy for values of the Reynolds number less than 4. 
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SUMMARY 


This paper is the first of two dealing thoroughly with the presence of an airfoil in 
shear flow he method to be adopted is illustrated by means of a simple example 
which also serves as a comparison with the results of other authors. The case of the 
general cylinder fixed at an arbitrary point in the field of flow is considered next, 
and the components of the force and the couple are obtained. The transformation 
s then modified so that the airfoil is either thin and curved, or else tapering and 
straight, and the general results reduced accordingly. 


1. THE problem of shear flow was first considered in some detail by 
Tsien (1), who showed how the complete stream function could be obtained 
for combinations of various types of potential streams with shear flow. 
The forces and couple on an airfoil in this type of flow were obtained by 
Kuo (2) from the modified Blasius formulae, but, in general, I have 
evaluated the forces directly from the pressure equation, derived from the 
equation of motion of the fluid. In a more recent paper, Coombs (3) 
evaluates the complete stream function for a particular example of an 
airfoil in shear flow by adding certain terms to the undisturbed stream 
function, which then satisfies the necessary conditions on the boundary 
of the airfoil and at infinity. This method, while applicable in the case of 
a stationary airfoil, does not so easily admit of a solution of the problem 
when the airfoil is moving. 

The fundamental condition to be satisfied by a stream function is that 
it shall always satisfy the Helmholtz equation, which for two-dimensional 
motion reduces to Dw/ Dr = 0, w being the vorticity vector and 7 the time. 
Thus, defining the complete stream function as ¢ = y%)+4,, wp being the 
undisturbed stream function, we see that the above condition will deter- 
mine the form of y,, so that adding several of these terms to %p, it would be 
possible to evaluate the relevant constants from the boundary conditions. 
Such is the method of Coombs. However, an alternative method making 
use of the Poisson integral is very suitable for the moving airfoil problem 
and does, in fact, give functions of the desired type in the case of shear 
flow. The functions found by this method are such that they reduce to 
specified value on a certain unit circle. If then we adopt a transformation 

z(t), which maps the boundary of the unit circle |f 1 in the ¢-plane 
on to the boundary of the given airfoil in the z-plane, on which & has a 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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specified value, we can apply the above method, the functions so obtained 


being such that, when added to the original undisturbed stream function, 
they satisfy the boundary conditions for the stream disturbed by the 
presence of the airfoil. The forces and couple on the airfoil may then be 
evaluated by determination of the pressure distribution over the surface 
of the cylinder and use of the Blasius formulae. 


2. The complex function Q 
In (1) Tsien deals with the problem of the stationary cylinder in shear 
flow by using the real stream function %, defined in the usual way in terms 
of the velocity components u, v, by the equations 
Cus Ors 
“= ——, v : : 
cy Ox 
which satisfy the equation of continuity in two dimensions, 
ou , ov 
ben oe @, 
ox oy 
While the method he uses is applicable to the case of stationary bodies in 
the stream, it is advantageous to use a complex function method when 
dealing with the problem of the motion of a cylinder through such a stream. 
To illustrate the method in its simplest form, I will first consider the 
stationary cylinder in the stream and show that it yields results in a closed 





form for any shaped cylinder, provided the region outside the cylinder may 
be transformed to the inside of a circle. The results for the circular cylinder, 
which take a very simple form, are then shown to agree with those given 
by Tsien, and the results for the Joukowski airfoil are given in closed form 
throughout. 

Following Tsien, it is assumed that the undisturbed stream function 
before the introduction of the cylinder is %, where V7, = c, a constant 
for shear flow. Since the fluid is assumed to be non-viscous and incom- 
pressible, the vorticity W associated with the flow must either remain 
constant or be such that it satisfies the Helmholtz equation, which in two 
dimensions is 


Dw 
— Q), 
Dr 
or, using the stream function %, 
ob 0 Ob Ole, 
ee V-s = 0. 
CX CY CY CX 





Putting / = J +4, and using Vs, = c, the condition becomes 


E (us, | yb) ei (us, ; yh) “|v, 0. 
Cx cy CX | 
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cy 
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The simplest solution of this equation is, of course, V*s, = const., which 
means that the insertion of the cylinder introduces extra vorticity into the 
fluid. But the effect of the cylinder is imperceptible at infinity, and hence 
the increase in vorticity is contrary to hypothesis, and so the simplest 
solution of (2.1) is V°s, = 0. Using this result it is possible to find a 
function & = y.+-y, satisfying all the boundary conditions of the problem, 
and also to define a complex function Q of which ¢ is the imaginary part. 
We may note here that since the flow is not a potential flow, the corre- 
sponding real part of Q is not a velocity potential. 

The method used to evaluate Q,, given ys, and the boundary conditions, 
is one which makes use of the Poisson integral to express the value of an 
analytic function in a region in terms of an integral round the boundary 
of the region, when the real or the imaginary part of the function is given 
on the boundary. This integral in its complex form has been effectively 
used by Sokolnikoff (4) in connexion with elastic problems. For a region 
inside a circle he shows that 

F(t) . J (9) 


do, (2.2) 


where F(f) is an analytic function of t, which is such that its real part has 
the value /(@) on the boundary of the unit circle y, ¢t and o being respec- 


tively points inside and on the circumference of this circle. 


3. The complex function 0, 

In applications of the above method we shall have to find a function 
given its imaginary part on the boundary of a unit circle. This case can 
obviously be dealt with since the conditions reduce the problem to that 
of finding a function —7: F(t) which is such that its real part is defined on 
the boundary. We also wish to apply the method to the region outside 


a cylinder of any shape. To do this, the region outside such a cylinder is 


transformed into the region inside the circle |t | by means of the 
transformation z = z(t). The boundary condition then being expressed in 
terms of ¢ only, the above general solution will apply, since # = 1/t on the 


boundary 
In the particular case of the stationary cylinder in the stream, the 
boundary conditions take a simple form. Tiey are: 


(2) that wu must reduce to ub, 


, at infinity where the presence of the 
cylinder is imperceptible, and 
(b) that y% must reduce to zero or a constant on the boundary of the 


cylinder, which must be a stream line. 
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Writing ¢ = g)+y,, the latter condition reduces to 


re(iQ,) = wp. (3.1) 
Using the transformation z z(t) to transform the region outside the 
cylinder to the space inside the circle |f 1, we have that 2 z(t), and 
that on the circle¢ = 1/# = o. Thus, on the boundary, , may be expressed 


in terms of o only, and the condition (3.1) is therefore sufficient to deter- 
mine 7Q, by the integral method. It yields, moreover, an analytic function 
of ¢ at all points within the circle, which is equivalent to an analytic 
function of z at all points in the region outside the cylinder in the z-plane. 
This means that the condition at infinity, namely, that %, should vanish 
there, is necessarily satisfied, and that only positive powers of ¢ will be 
given by the integral. Thus all the conditions of the problem are satisfied, 


and ier do 
102, (t) : {7Q,(0, ¢)—71Q, (0, &)} ‘ (3.2) 


4. The forces and couple on the cylinder in terms of the complex 
function 0 
The components of the force, X, Y, on the cylinder, and the moment 
of the couple [ about the origin, are given by the integrals 
Y -+2X | p dz, r re | pz dz, (4.1) 


( Cc 
taken round the boundary curve C of the cylinder, where p is the pressure 
at any point of the surface. 

Returning to the equation of motion of the fluid, we have that the 
pressure distribution throughout the fluid is given by 


P L(u?+-v?)— Cy, 
p 
where the constant C is determined by V*xs, = C. 

However, for a stationary cylinder, 4% is zero or a constant on the 
boundary, and hence gives no contribution to the forces and couple. The 
forces are thus given by 

Czy 4 - kp | (u?-+-v?) dz, 
es 
which in terms of 2 and ¢ can be written as 


Y+ix =" [ eQ ef) AQ A 
y)la alla at 








dz\-1 _. 
dt. 4.2 
("1 7 


the double change in sign in transferring to the circle y in the ¢-plane being 


necessary since integration round the boundary C of the cylinder in the 
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--plane in the positive direction means integration round the circle in the 
-plane in the negative sense. 


Similarly, the couple is given by 


"TQ afeQ A] fd2\_. .. 
r = re (—4p)} |< il ( z(t) di. (4.3) 


5. 2 for a circular cylinder 
In order to facilitate immediate comparison of results with those of 
Tsien, we will first consider the horizontal type of shear flow that he uses, 


| umely ) 
u ly kK? : v 0, 


| . 
| ( 


measured with respect to a fixed set of axes in space, Uy, K, ¢ being con- 
tants. Considering now a circular cylinder of radius a fixed at the origin 
ff coordinates in the stream, whose axis is normal to the field of flow, we 

iy transform the space external to the cylinder into the interior of the 


nit circle |f | by the transformation z a/t. Thus we have that 


faif/l 1 a*kK (1 1\*) 
ub =|-+ =I }, 
a oe ( i} 8c (; t} J 


so that 42, is given by 


nd hence by the integral method, (3.2), 


iQ,(t) = Uy| —ait 4 in| 
| 4c J 


so that the complete function Q is given as 


To compare this result with that of Tsien, it will be necessary to write 


c/2r)e whence, evaluating uw, we have 


c* . K fe* . é | 
l AG }sin 6 sin?é 4 - COs 26 | 
i 2\¢ 32r? 
rhe difference in sign between the two results can be accounted for by the 
ict that Tsien defines his stream function ys as being such that u = oyb/ey; 


snowing that. in fact, his us corresponds to our us. 
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6. The forces and couple on the cylinder T 
The components of the force on the cylinder are given by give 
re l ion l . Ka? :.4 
yeix 28 f vgle(s—)—ik@(_ 1,1), 
2a} | t? Ze\ @'t J 
; , , Kae ; l = 
la(l t?) +2 . ‘aa a }j ae whe 
” t}) inte 
. i. ee en 
which reduces to Y+iX 2rpU3 K —, vot 
C ars 
pars 
. , rg pr a" 
whence X 0, ) 27pU2 K 
r 
which displays a lifting force only. This again agrees with the results of x 


Tsien when we put a = c/2, the radius of his cylinder. 
Similarly the moment of the couple is given as the real part of the 


integral so 
e Ka? 9 r 9 Q) fe 
, vel; Ka 1 a ,Ka?l ;Ka* | 
2P oe a aS : : and 
3 | 2¢ 8 # 2c ft 2c | 
y (a 
Ka? 1 . Ka? . Ka? ,,) dt : 
— L@-4 t—at?—1 | i (D 
zo. t 2c 2c ft (c 
which is zero. : 
The 
7. A stationary airfoil in shear flow the 
Having established the method of procedure for the case of stationary 
ee : : ; iQ 
airfoils, we now introduce the shear flow with which we shall be concerned — “1 
in the remainder of this paper, namely that shear flow whose vertical 
velocity is a linear function of the distance along the horizontal axis: 
u 0, v Uji+ke], ; 
whic 


where k is a constant, measured with respect to axes fixed in space. There 
is no loss in generality in using this particular type of flow, since the  g ¢& 
' ; ; oe ois 
results for the stationary airfoil in the shear flow defined by 
u = Uj1+ky], v= 0, 
may be deduced very readily from the general results. 
With respect to any set of parallel axes with origin at point (%,¥) | Thu 
referred to the original axes, the above flow may be written as 
u = 0, v = Uf A+ke’], 


where A = 1+k%,, and 2’ is the new coordinate. The point (%p, J) is 


determined by the position of the cylinder in the field of flow, and is in 
fact the centre of its normal section. In this way it is possible to measure 
the forces on the body at all points in the field, and to plot their variation 
as the body is placed at desired points, though the airfoil must be held 
stationary at all positions. 
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The two-dimensional airfoil to be considered is the general cylinder 


given by the Joukowski transformation 


when the space external to the airfoil in the z-plane is mapped on to the 


dt 
interior of the unit circle /|f |. If the cylinder is inclined at an angle 6 
to the horizontal, then referred to axes O'X’Y’ through its centre and 
parallel to those fixed in space, the transformation becomes 
a c* = 
2074 7 et”, (7.1) 
t 2 +a/t 
ults Now the original stream referred to the O’X'Y’ axes becomes 
u—0 v = U{[A+kz’], 
tu , was ; as = A 
so that uy = U,| Axv’+4kev’?], which is the imaginary part of Q). Defining 
Q for the complete disturbed flow as Q = Q,+Q,, the conditions for Q 
) 
d { , are 
a) Q Q, at infinity, 


b) & is a constant on the boundary of the cylinder, 
c) VQ, 0 at all points. 
The second condition leads to re(iQ,) = U,[ Ax’ + $ka’?] whence, employing 


the integ1 ul method, we have 


jonar’ u ff { | a C“o_ | Cc? 
ETT S2; | ) Te Jet 9 pede et) + 
oi 7U . | 2 Co al Lo 29 7- ao 











ert 
i a co c2 . : d 
’ {zo et? Z9 + ao + eto J , 
8 ao a+t+Z oa Z9+ao Jo—t 
Ther which reduces to 
2 . 472 
ce the Q(t) ai IA el P ate-19| 11 h c't e2i0 1 
\2 ja+<zot 8 | (a+z,t)? 
st ; i cZat act? ac4 | 
(a-t- 20Z,t)e-*" 2{ Zot id | ———— : 
A+zt a+zyt  (a+zt)(a®@—Z29Z)) |) : 
Thus the complete Q for the disturbed flow is 
, / 1 
Q) iV | (z } ate “| | ' 
| t 
rf : | | = 4 fs 4%, cv ei} __ 
d is Sil" t | atzot "2 e+ 
easu k c4t2e2 242 az. t PY.) 
ort (a*t? + 2azZ, t)e-?" 
riatl t 1 (a 0 ) 
ye he 949 
a yf act” | ac} | 


+ 2t)(a*- 
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8. The forces on the cylinder 


Using the formula (4.2), we may evaluate the relevant terms from 


(7.2) and obtain the components of the force on the airfoil as 
Y’1ix’ = — 2 U8, iO 
rT 

| (> a gn, back , back? kare | 

© \n<1 n 2(at- Zu) 2(a Zt) 2(a- zyt)(at =, )(a2 9 25)| 

<| > B,t"-* katt L Broa eg katct ) 

\n&1 2(a-+-2ot) 2*(at+-%) 2(a+-29t)(at+-Z)(a°—2yZ,)} 
(a+ Zt)? di 


| [oc (a T Zyt)?| 


where the coefficients are given by 


ka? ka ha 
. 20 Q.. i0 te eC ed P. 
Xy 5) é P53 X9 Aae gee <of >a 9 “0 P43 
ka? ha ha 
. Q. Mipcom 3 » ptt O. 
X3 > P33 X4 Aae™ + > *at> 70e" Po; 
ka? , 
© 2i0 *] , 
X5 = —>-€ Py. (8.1) 
The integral being taken round the unit circle /f 1, the only singu- 
larities giving non-vanishing residues are the poles at ¢ = 0 and ¢ 2,0. 


Evaluating the integral term by term, and writing p? = a2—z,2,, we 
have, after some reduction, 


pan) Ses —_—" c* .. CO" 
Y’+iX’ = wipU2k Aal 2a4° e204 © 9-210) 4 

| a a 
hk: ‘ c ; oy) 
~ 1 9g 2s Ct «. ri. care 208 1s. 20 * 01s 6 
5 | maze ze) +6 (Ze +2 e-7) — (age +2ge-")| 1s 

p? 

which, being entirely imaginary, gives only a forward force X’ as 


-— 2apU2 kl Ala? +-c? cos 20] + 


\ 


4 “4 
mr A | 2a2+-¢2_S ~|a) cos 6 — | 2a2—c? — e Y,)sin || (8.2) 
2 p- p- 


These results are not entirely unexpected, since the presence of a force X’ 
only indicates a force urging the body across the stream. We can draw an 
analogy between the above and the case of a cylinder in an ordinary 
potential flow. If in the latter the airfoil is stationary, no force is ex- 
perienced by it unless a circulation is introduced round the cylinder. In 
this case, for a positive circulation, there will be a force on the airfoil 
tending to urge it across the stream, the direction of the force bearing the 
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same relation to the direction of the stream as does the force X’ above to the 
fror direction of the shear flow. Thus we see that shear flow in the field of non- 
potential flow is analogous to a circulation round the cylinder in the field 
of potential flow, which, since circulation and vorticity are closely linked 
when dealing with airfoils in potential flow, agrees with our knowledge that 


shear flow is a continuous distribution of vorticity throughout the fluid. 


20) 9, The couple on the cylinder 


Using the equation (4.3), the moment of the couple on the airfoil about 


mul the origin at its centre is given as the real part of the integral 
0 a 
. Ee dias Lar2 kac2t ha2c4t? 
les 2(at+Z,)  2a+tzt) 2a+z9t)(at+zZ,)(a*—z,_Z,)) 
= ast zc’ 
($9 pn kac*t hac . ka .. 
\n=1 2(a+z jt) 2#(at+2Z,) 2(a+2,)t)(at+zZ,)(a?—Zz,Z,)) 
(a t) |< 2 (a Nf 2] 
a —> at, 
ae f 72 (a+29t)?] 
Q 1) ] | 


where the values of «, and 8. are given by (8.1). 











singu This integral, on evaluation, becomes 
— ee 2] 14 222181 9 4h] 2q2e, et 4 c2z, ci8 4 ¢4~9 e-i8 
Cn. WE » f 0 0 0 9 
) | p? 
6 2 2,6 
» 9 } 9.9.9 ( c“ ac 
}:2| @222 e2 040418 1 9q2e2e2i8 52 — c422 
0 a9 252 4 2 ¥ 
p (c“29— p*) \ p 
1 6 
( im( 
my 2 oe, Fs 452[ 2 Py | 
|) => j | i 0 a ti | Zo] za « il) | 
f | rm pP i pp 
The real part of this expression gives the moment of the couple as 
mpl 0 
4 4 
. . ° ° ( 6 ( 
12 422 sin 204 Ah rosin 6{ 2a ce -| + y, cos 0{ 2a2—c? — — 
| p- p- 
Sve J.2 | (6 
5 > =. ° o ,@ 2 ¢ . 0 
= | A*| | eo— Yo [SIN 20 +- 229 Yo COS 26) -ct sin 40+ 2a%c? sin 20+ 2— x9 Yor 
| , 
. = | 
oree A ' , ,; 
DatalOna p65 DP p2l 2 o 
raw al aC XH Yo =O°XG Yol = [ae yi | p”) 
ps Ise 1] 242 4 [ n252 | 2642 4 
rdinal P| O'#9 P \Lo <0 Pp Lo “0 P Il “0 r | 
: Y72710 2 Qe = ; 14,2 52 8 97105 442 52 8 
IS eX | = a 220 29 )Xo Yo ‘ “0-0 pP ) ~( Xo Yo = be p ) _ 
j 252 412[ p2,2 42 4121 p22 4)2 
ler. hh pr “0 P ak “0 pP | P |< “0 Pp | 


40125 25 Xq Yo(C*|x>—Yo|—P*) | | 
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[t is usual to take z) small so that its squares and higher powers can be 
neglected. In this case we obtain for the couple 


4 mpU? . 
- 4 4 
12 42¢2 sin 20+ Ak X_ Sin 6| 2a?+-c? — —]} + y, cos 6| 2a2—c? — 
| a= ‘ a2 
k? P oe 
t= [c# sin 40@+-2a"c? sin 26}. (9.2) 


We observe that, in this instance, for 6 zero we still get a small negative 
couple on the airfoil. 


10. The effect of thickness and camber 

Having obtained the forces and couple on the general cylinder, we can 
now examine these results in the special cases when either the thickness 
or the camber of the airfoil predominates. This is effected in the usual way 
by making assumptions as to the nature of z,, the centre of the circle of 
transformation in the z-plane. 

Thickness. We may examine the case of a straight tapering airfoil by 
assuming that z, is real, that is y, = 0. In this case the horizontal force 
on the cylinder becomes 

rr ) a) { 9 9 ) i: Yy2 9 | 
X 27 pU?2 k\ A[a?-+-c? cos 20] + — | 2a2+-c? — x, cos 6}. 
i » € on 
| 2 a*—X5 
For an ellipse, 2, is zero, so that the increase in the horizontal force caused 
by the difference in the thickness at the ends is 


4 

D2 2 ( 

a To 9 9 
a*—25 


mpl f k? XV COS 6. 








Now a >c > a, s0 that the above term changes sign with 2, and cos@, 
but for a given value of 2, it maintains a constant sign for values of 6 


given by —}$z7 < 6 < $n, there being a maximum of 
= ct 
rpU? k*| 2a*-+-c? ______|a, 
a?— 72 
at 6 = 0. Since the factor A 1+-k%, is absent, then the extra term is 


not a function of position, and is constant for all points in the field of flow, 
though the complete force X’ is seen to be a linear function of #p. 
Similarly, from (9.1), the moment of the couple becomes 


r 

‘i re 6.08 " ‘ 9 ‘ € 

[ mpl 2/94 22sin 20+- Aka, sin a( 202 +-c? — — ) 4. 
a?— x? 


\ 


2 
~ 


9 


ke? re i ee 
a (a2? sin 20+-c4 sin 40+ 2a?c? sin 26), 
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and the extra terms introduced by the tapering of the airfoil are 


r | az—x2 


npU2) Aka | 2a2-+-c2 





k?a?2? cos 6\sin 0. 


| It is evident that for lar < 6 ir, this term changes sign with sin 6; 
and since 2’y is very small, the second term is almost negligible, so that the 


l>, and are 


° 


maxima occur at inclinations @ 


9.2 
1 
( 
amol’= Ak| 2a*-+-c? =|Xo 
tive “ss 
The presence of the function A in these extra terms indicates that they 
depend not only on the shape of the airfoil but also on the position of the 
eylinder in the field of flow 
eas Camber. To examine the forces on a thin curved airfoil, we let z, be 
ness entirely imaginary, that is, 2) — 0. In this case the horizontal force 
Wa becomes 
sis X 2rpl 2 I) 4 a*+-c? cos 26] + — 2a*-+-c* 4. —___.. yosin 0\, 
| 2 a?— ye 
lL by and the additional term introduced by the camber of the airfoil is 
ores oA ; 
mp Us k*| —2a*+-c* + —— | yp sin @ 
a°— Yo 
which is generally negative, but becomes positive for negative values of 6. 
: In spite of this, the horizontal force on the cylinder is positive, and is a 
used ss . . ae 
inear function of Z,, the distance of the centre of the airfoil from the 
fixed Y-axis 
Similarly the couple on the airfoil becomes 
o 
08 8 7pU?2!2A%c? sin 20+ Aky, cos 0} 2a?—c? : si-4 
| a Yo 
ot 
| —a?y? sin 20-+-c* sin 46+ 2a7c? sin 26}. 
ind the extra term due to the curvature is 
aia sil o —— 
sella mpU?2 k\ A | 2a2—c? | —ka*y, sin 6)y, cos 8, 
flow a°— Yo 
i linear function of @. This is a negative quantity which increases the 
negative moment of the couple on the airfoil. 
The results for the case of the stationary airfoil in shear flow whose 
velocity varies linearly with the vertical height may be deduced from the 
) ibove by putting @ Lar . where @’ measures the inclination of the 


airfoil to the new X-axis. In this case, calling the components of the force 


it Ee 














$18 TWO-DIMENSIONAL AIRFOILS IN SHEAR FLOW 





and the couple on the cylinder X”, Y”, and I”, we would have X” — Y’, 
; =J2",eua i” Tr; also A” 1-+-kyo, so that for a general cylinder 


in this type of shear flow, we have the following results. For the forces, 


Y” — 2npU2k> 


” > > , k > > a . , 9 9 c 
I 4 [a*—c? cos 26’|-4 2a?+-c? x, sin @ 2a2—c? y, cos 6’ || 
= 5]“0 3]0 
p- p- | 


and the moment of the couple for small Zp, 


"7 ro "oO s , , , > ? ce 
I mpl 32 2c? sin 20’+-A t| x9 0080 (20° ¢*.. ; 
az 


9 


a* 


yosin 0 (20° e "| |e [et sin 40" 2ate? sin 26"), 

Although the results obtained above are of interest in that the example 
considered, the general cylinder, may be reduced to various special cases, 
the paper itself is intended mainly as an introduction to a subsequent 
paper which will deal with moving airfoils. On examining the integral 
method used here it is abundantly clear that the functions so obtained do 
comply with all the conditions of the problem, and are in fact determined 


by the conditions on the surface of the cylinder, so that in the case of 


moving airfoils it will still be possible to determine the complete complex 
function Q, since the value of % on the boundary will be known. The 
problem of the generally moving cylinder with arbitrary circulation in 
shear flow will be solved in the second of these papers. 

[t must be emphasized that the solution of the hydrodynamical problem 
results from only one of the several mathematical solutions of the equation 
(2.1), but it is necessarily unique on the fundamental assumption that the 
presence of the cylinder causes no variation in the total vorticity of the 
fluid. This solution will therefore be used in the subsequent paper on 
moving airfoils. 
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“— JETS PAST FLAT-PLATE WINGS. II 
P By M. HOLT (Department of Mathematics, The University, Sheffield) 
veceived 17 April 1951] 
SUMMARY 
The method of determining the pressure distribution on delta wings in a non- 
iniform supersonic stream developed in (1) is extended to include wings on which 
26')\. subsonic conditions apply ng part or the whole of the leading edge. 
| 
imple 1. Introduction 
ile In the first part of this paper we used the linearized method of conical 
quent fields to find the disturbance in the stream from two adjacent supersonic 
tego) jets due to a delta wing with its apex on the jet boundary and with leading 
ed do edge lying entirely outside the two Mach semi-cones. We now find the 
nined disturbance when the leading edge in either jet, or in both jets, lies inside 
sse of the corresponding Mach semi-cone. Three cases have to be considered: 
mplex (a) leading edge inside the Mach semi-cone in the port jet, outside the 
The Mach semi-cone in the starboard jet: 
ion il b) leading edge outside the Mach semi-cone in the port jet, inside the 
Mach semi-cone in the starboard jet; 
oble c) leading edge inside the Mach semi-cone in both jets. 
iatior In each case a procedure similar to that set out in (1)is adopted. Firstly, 
at the the flow outside the Mach semi-cones is determined to give the boundary 
of the conditions on the Mach semi-cones themselves. The flow inside the Mach 
er O1 semi-cones is then related to one unknown function which must satisfy a 
certain integral equation. Wherever possible the notation of (1) will be 
used here. 
2. Wing with leading edge inside the Mach semi-cone on the 
port side 
The plan of the wing corresponding to case (a) is shown in Fig. 1. The 
p. I origin of rectangular Cartesian coordinates is taken at the apex of the 


wing with the z-axis in the direction of the main stream and the 2-axis in 
the transverse direction. The velocity, density, and Mach angle in the 
main stream are W, p, » respectively in the port jet and W’, p’, w’ in the 
starboard jet. AOA’ is the leading edge and OB, OB’ are the traces of 
the Mach semi-cones. The angles of sweepback are denoted by f and f’ 


) 


where in this case, 8 t7—p and pf’ < }$a—p’. 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 











‘is the wing leading edge, 


of the vortex sheet, ¢ }’ are the traces of the 


Mach semi-cones. 
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In the port jet (a 0) we employ non-dimensional variables x,, y,, 7, 8 


lefine d by 


tan 2) xy r cos 6, 


tan j2) Wy r sin é. 
Corresponding variables in the starboard jet, x;, yj, 7’, 9 are defined by 
similar transformations with « replaced by pv’. Conditions in the (7, #)- and 
planes are shown in Fig. 2. The semi-circles 4 BA’, DED’ are the 


traces of the Mach semi-cones and LCF is the trace of the wing. 


. . 1 F . 
Conditions outside the Mach semi-cones 


\s in (1), only the disturbance in y > 0 needs to be considered. The 


solution outside the starboard Mach semi-cone and the boundary condi- 
tions in the (7’, @)-plane are exactly the same as in (1), section 3. On EF 
id within region 5, v W's. Within region 6 the disturbance is 


related to that on the section AD of the vortex sheet. Elsewhere the 
disturbance is zero. 


In the port jet the disturbance outside the Mach semi-cone is zero. 


7 ] ~ 
(Conditions side Tire Mac /, S€7NL-CONES 


Inside the Mach semi-cones the disturbance satisfies the boundary con- 
ditions already found on the Mach semi-cones themselves and the addi- 


tional conditions, 


W'n on CER, 


Wx onCl, 
Ww 0 on BL (giving continuous pressure across the vortex sheet BL). 
The disturbance is defined by analytic functions G,(t), @,(t), G5(t) and 


Gt’), Gv’), G3(t’), in the port and starboard jets respectively where 


lt sin 6 {7(1—r?)? cos 6}/r, (1) 
Lit K sin 6/r’ +- KSi(1— r?)! cos 63/7", (2) 
ind 
IG, /dt:dG,/dt:dG,/dt (1—#?)! i:ittan p, (3) 
1G, dt’ :dG dt’ :dG@y/dt’ (l—K#?)?: —i: it’ tan p, (4) 


where principal values are taken for the square roots. 
\ll these functions involve the unknown function defining the distur- 
bance on the vortex sheet CDA. To determine this we set up an integral 
juation expressing the condition that pressure and normal velocity are 


continuous across the vortex sheet. 
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Solutions in the t’-plane 
The solutions defining the disturbance in the ¢t’-plane are given in (1). f 


equations (16) and (17). They are 


1 
dG, ci [ 2sP\(s)ds | 2W’atPitanp (5 
dt ot | {2 —s? 7(t2—t"2) ) 
0 
1 , , » P , 
dG c(1— K2t’?)! } 2sP1(s)\ds 2W ‘at? (1 Kt") (6) 
dt’ mt’ tan pu t'2— 8? at’ (t—t'?) 


yf 


where u’, v’, w’ are the real parts of G{(t’), Git’), G3(t’) respectively and 














Si 
w’ cP(7’) cPi(r), 0< r(=r) <1 = whentl’ r—QO. 
Solutions in the t-plane 
Fig. 3 shows the transformation to the ¢-plane from the (7, @)-plane of T 
B 
Fy 
07% || L 
WI 
fin 
Th 
B 
7 A > 
> 
0 1 E 
oi oe wh 
Fic. 3 
Fig. 2. Ifr, = cot wcotf is the value of r at L in Fig. 2, the corresponding 
value of ¢ is t= 17, where 7, r,/(1—r?)!. | 
The boundary conditions in the ¢-plane are 
w P(r) <= re 
t ra Oi. fal 
u,v,w = 0 s<r<-e rh 
\ 
v Wx t} <7, on the imaginary axis, 


- 7, on the imaginary axis. 
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The required solutions are most easily found by the use of a step function 


1 for P(r). Accordingly we find the solutions satisfying 
| Q s | . 
t r+On2, 
4.2 () S x } 
7 Wr f Ty | - ° . 
on the upper imaginary axis. 
{ 8] Ty 
Thus. on the real axis 
(6 
re dG, dt 0 t - T1> 
im dG,/dt = 0 |t| > 7. 
and 
Since w is finite at infinity IG, dt O(1/t*) for large t. 
The function satisfying these conditions is of the form 
dG. Ai Bi (7) 
i ; 7 
dt t*-+-79)#(t?—8*) © (+73) 
une Ol To find A we integrate this relation and use the condition 
re G:(t) 1 (0 t =< #@). 
From (4) and (7) 
dG, Ai Bi (8) 
dt ttan p(t? T?) (¢2— 8) ttan p(t? + Ti)? 
which, when integrated with use of re G(t) Wx, |t} <7, gives B. We 
find > > 
: A 28(77-+- 8°)? /77, 
B 2Wary tan p/a—2rH{ (77-4 g?)} 7}/78. 
The solutions corresponding to the full boundary conditions are 
1G, i(1—t#*) dG (9) 
dt ftan p dt 
where 
1 
dG. ) * 28(7?-++- 87)? P1(s) ds 
ondu Vv , F _2 | 2 2 I 
dt - Tt)? f s*) 
Din? | 7 +s?)§_7,'P'(s)ds  2iWar3 tan p (10) 
a(t Ts : Ss a(t? + 72)! 


The inteq al equation JO) P 
We now substitute equations (9) and (6) in the relation 
dG dG, t ) 


O71). 
re . re . ‘ for 0 < r < il, 
dt if t r—Oi } 
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which expresses the condition that the normal velocity is continuous across a 


the vortex sheet. The following integral equation for P(7) is then obtained: — | p 


1 
(1—#)! [ 2s P1(8)(7}+-s°)' ds 
) 2 


re : " = > — 
tan p(t°+-77)? . t s2 
0 
) l » 
275(1—?#*) [ P1(s){(7§-+-s?)? —7,} ds 277 Wa ) . 
tan p(t?-+-77) S (t?-+-77) 


1 
c(1— K¢?)! [ 2sP\(s)ds | 2W’at?(1— K?)! 


= — Q (11) 
t- 8° Bey 


tan pe 
0 


(t r+0Oz, 0 = 2. 


We may define P'(7) for negative arguments as we please. If we define ’ 
it as an odd function we may replace (11) by 
( 
reH(t)=0 ((=r4+0, O<r<c)), U 
where 6 
1—?? > Pl(s)(7?2-+- 8)! ds 
H(t) ( : ) | (S)( 1  ¢ 
tan p(i--T7)*? . s—t 
> 1 € r 
7*(1—??) | P1(s){(7?-+-s?)!—7,\ ds , 27? Wa(1—#*)! 
tan p(t?-+-77) S (t?+-72) 
1 
r 1 1 ry , r 
c(l1—Kt?)! ¢ P\(s)ds | 2W’at?(1— Ke) (12) 
tan pe J s—t te % ; F 
a 
We make H(t) single valued by cutting the ¢-plane along the real axis 
between —1 and +-1 and along the imaginary axis between — iz, and (7. 
Then re H(t) 0 for all real ¢ and H(#) is regular in the cut plane except 
at t -iz,, where it has singularities of the form 
a(t--ir,)-*+-b(t+i7,)? , 
and at ¢ -#,, where it has simple poles. At infinity H(t) = O(1/1). We 
therefore take a 


Ait ; Bit Cit (13) 


ae an Oe) 
(t2—t?) © (#+-72) (tg —t*)(t*-+-77) il 


where A, B, and C are real constants to be determined later. Equation 


(13) defines H(t) in the whole of the upper half-plane and can be continued 
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cross ynalytically to define H(t) in the lower half-plane. Then taking imaginary 
ined parts when f rt+0O2z, 0 1 kK, we obtain 
pI tan | B C 
r=) c( 1 K?;?) lt? 2 (r= Ti) (t? r?)(r2 7?) | 
(14) 
nd when rtOi. 1/K l 
4 , 1 re , r 1 
7(1—r?)! Pl(r)__ e( K2r?— 1) 2 * P's) ds__ 2W'at?(K?r?— 1)! 
tan tan P s—r t? r2 
. 1 
1] 
1) Bi C) - 
—" a oe reer (15) 
t2—y2 ' (y2+-72)8 ' (42 —r?)(r2+-72)! 


where 7 denotes the principal value of the integral. 

Equation (15) is an integral equation of a general type considered by 
Carleman (2) and may be solved by extending a technique employed by 
Ursell (3) in solving a simpler equation of the same type with constant 
oefticients 


Considel 


P\(s) ds 


((t) (i—J 1 —?¢”)!+-c(1— K#?)!! 
« 6 t 
1 
Then, comparing with (14) and (15) we see that for ¢ lrtO,0<r<I1/K 
A B C 
m Q(t (a tan php} —,—, + = aT a 
) tty y (7°-++77) (t5" r-)\(re- T?)}} 
ind for 1/A i 4 r+or 
| ax 
' ' ( A b ( 
id ir im V(t r re) tan i =; =_-] = = : = cee . 
| a iy Ti) (t-—r-)(r- T{) | 
2W ol? tan p( Ar? 1)3(1—7?) 
sy r2 
For f yO ‘ 
im O(f QO; 
W, n Y(t) 
ulso Q(t) = O(1/#) at oo. 
9 Thus Q(t) is defined everywhere in the upper half-plane by Poisson’s 
i integral . 
1atiol | Ort) 1m Qu) du 


7 u—t 


inued . 
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Hence, 
1 
| P1(s) ds (1—#*)! 
; s—t a{(1—??)?+-c(1— K #7)! 
1 
1 
r utann {( A B ; C \y 
; a oo + au — 
| (1—w?)*(w—t) ((t2—u?) © (w?+-72)! © (#2 —u?)(u?+-7?)}] 
1 
VK 1 5 : : “ 
| } 2W ot? tan p( K2u?—1)! du 16 
} oe (1—w?)? (#2 —u?)(u—t) ee 
1 1/K ‘i 
and when 1/K <r<l, 
nP\(r) a : yr) = oocens oh, . A hae B = 
(1—r?)—c?(1— Kr?) | (L—r?)! (t2—1? © (r?+-72)! 
c )  2W’at? tan p(K?r?—1)#]  e(1—r?)!(K?r?2— 1)! 
(#2 —r?)(r2+-72)!] (1—r?)!(t2—r?) ' (l—r?)—c?(1— Kr?) 
1 
L|F | —— [_ A | B mi C \ de E 
J (L—u*)(u—r) t?—u? © (u? +79)! (t2—u?)(u?+-79)!) 
1 
VA 1 


7 [ ; [ 2W'at? tan p(K2u2—1)! dul (17) 
I (1 u?)i(u r)(t? u?) | . ‘ 


1 1K 
The coefficients are determined by the condition that the singularities 
on both sides of equation (13) must cancel out. We find 


A+(C/(t2+-72)! 2W’ at,(K2t2—1)}, (18) 
1 pn 
, . ] 7 s)§(72-+ s?)i— +, ds 
B 7(1-+-72)!| 2Wa— | ci a a (19) 
tanp , 8 
-1 
> Y 9 i 1 9 9\1 
3B Cr, (3—74) Wy — l P1(s){(7?-+-s?) vi} ds] ; 
87,  (t?-+77) 8(1-+-7%)! tanp , Ss | 
4 : 
72)! S 1 S Ss 
(1+-77) sP1(s)d (20) 
tanux J (s?+-7%)! 
1 


where P1(r) is given by (14) and (17). 


The pressure distribution on the wing surface 

When A, B, and C have been found, the unknown function in (5) is 
given by equation (16) and that in (10) is given by (14) and (12). The 
pressure on the wing can then be found by integrating (5) along the 
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imaginary ¢t’-axis and integrating (10) along the imaginary axis and noting 
the relations 
p pWw in the port jet, 
and p' po’ Ww’ in the starboard jet. 
3. Wing with leading edge inside the Mach semi-cone on the 


starboard side 


In case (b) conditions in the two jets are reversed. The solutions defining 


(16 the disturbance in the ¢-plane are given in (1), equations (14) and (15); they 
are 1 
dG; (1—#?) [ 2sP\(s)ds | 2Wati(1—#)! (21) 
dt mttan pe . {?— 8? '  mt(t2—t#?) 
0 























1 
dG. i f 2sP'(s)\ds 2Wat?itanpu 
and : | - (22) 
dt v {*—s* a(t5—t*) 
0 I/K 1 od 
f Ll 
u A Y; Cc D H —_ 
A E 
(17 
rities 7 To F 
r | 
Yo 
(19 / 
Fic. 4 
Fig. 4 shows the transformation from the (r’, @)-plane to the ¢’-plane in 
this case. 
At F,i it,, where 7, r,/K(1—r?)' and t cot »’ cot 8’. The 
boundary conditions in the t’-plane are 
(20 u cP(r) U0 ; 2 | | 
t r—Onr 
/ { 8) | / 
] Vi ‘a t Tr, | ° ° ° 
on the lower imaginary axis. 
u Q / To 
D) 18 The solut ions satisfy ing these are 
Th dG, i(i—K*t’?)'d@, 


the , , 
th | dt t’ tan p dt 
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where 


1 
dG Ci 28(72 gs”) Pl(s) ds 
dt a(t’? te) | (t’2 s?) : 


ier? { (721-52) 701 Pl(s) ds — 20W ar tan p 
2 5 UT2 2 i 2 I (24) 


Ss , a(t’?-+-7,2)! 
The integral equation obtained from the condition of continuous normal 
velocity across the vortex sheet may be written 


re H(t) QO (¢t rt+0:, O<r< 1), 


where 
1 
(1—K2)! ff PX(s)(7'2-+-8%) ds 
H(t) . = 2 | = i. 
tanp(l+7y)? , s—t 
1 
- 1 F ’ 317 4 
cre (1 K2#?) | Pl(s}{(r2 s?)! Te} ds ' 273 Wa I K*t*); 
tan p(t? 7) : S (t? te) ~<a 
1 
al " 
(1 {?) . P1(s) ds 2Was?(1 t*) (25) 
ee J 25 
tanp , —s - (fj t*) 
1 


Following a similar argument to that used in section 2 we find that 


A’tt . Bit Ct 


(@—#) ' (@+72)! ' (@—ey(e+72) 


H(t) (26) 
in the upper half-plane. Equating imaginary parts of (25) and (26) we 
obtain: 


For t rt+03n,0<r< IK, 


P\(r) rtan p { | B’ i. ( | 
mr{(1—r?)*+-e(1— Kr?) (B—r? © (rP+r?)t | (B—1?)(r?-+-72)!} 
7 
21) 
For t trtO, WK <r i. 
1 
. c( K2r2 Da Z [ P'(s)(7?-+-s?)! ds 
tan p(7? - T+)? s—r 


1 
er (Kr? 1)! [ P1(8){(72 |. 32)! Te) ds 
tan p(r* 
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We note that the integral occurring in the second term on the left of (28) is 


be denoted by D’, say. This second term may 


independent of 7 and may 
then be transferred to the right side of the equation which is now of the 
same form as equation (15). To solve it we use Poisson’s integral formula 


to define the function 


1 
(1—#?)!+-c(1— kK)! f P(s)(72-+-8?)! ds 
(0) 75 \ = ) (29) 
(|—i {-4-7.°) s—t 
1 
ind proceed as in section 2 to find P'(r) in terms of the four unknown 
constants A’, B’, C’, and D’. The condition that the singularities in 
equations (25) and (26) should cancel out yields three relations between 
these constants The fourth relation required is that 
P1(8){(72-+-s?)!—75! ds 
D ~s é (30) 
8 
When these constants have been determined the unknown function in 


24) is given by equation (29) 


The unknown function in (21) is then found from (31), (25), and (26). 


4. Wing with leading edge inside the Mach semi-cone in both 
jets 
In this case the solutions inside the Mach semi-cones, expressed in terms 
of the unknown disturbance function, are already known. In the ¢t-plane 


the solutions are given by equations (9) and (10); in the ¢’-plane they are 


given by equations (23) and (24). The corresponding integral equation is 
therefore 
re H(t) = 90 (t=r4+0, O<rcl), 
where 
F / s\(r2 s*) ds 
H(t “— 
7)" . roe 
l 
| | t | P! s){ (73 s?) 71} ds 
in p(t rs ’ s 
> 1 ’ 
273 Wa(1—?? c(1— K#?) * (72 +-8?)! Pl(s) ds 
Ts n p(t 7.2) p s—t 
1 


A P1(s) ds 
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Then H(t) is purely imaginary on the whole of the real axis, has singularities 
of the type a/f*+-b/@ at t +ir, and ¢ Liv, and H(t) = O(1/#?) at o, 
Accordingly, we take 
a ea Sn ck 
(f+)? (P4778)? (@+77)2(F +78)?  (+79)?(+72)! 
to define H(t) in the whole of the upper half-plane (cut along the line 
joining the singularities). H(t) can be continued into the lower half-plane 
by means of Schwartz’s reflection principle. Taking imaginary parts on 
the upper real axis we obtain: 

For 0< r:< 1/K, 


, ‘tan pu 
P(r) pee 
a{(1—r?)!+-c(1— K?2r?)} 
1 B Cc” D ) 
(r2-73)t  (r2rh yt (72-7) P)t OS (2-79) (7? +- 7) 
(33) 
and for l K < ? < :, 
1 
(1 ‘Tt . 22):) ll - @p [ (72-4 s°)! Pl(s) ds __ 
tan pu tan p(7"*-+7T.,")? 7 s—r 
1 
1 
cr(K*r?— 1) [ {(72-+-s?)—7,}P1(s)ds  273W'a(K?r?—1) 
tan p(7* t?) ‘| S (7? 73) 
1 
| B" ( D | 
} a oe ) > 19\3 9 i) 
(r?-77)P OO (r? rye) (9? ap)? re) (9? 77)* (7? +72") 
(34) 
To solve (34) we put 
1 
» tfoe 92) 7 'P1(e 9 
| (Ts ) ahi ( )d E” (35) 
4 » 


and then use Poisson’s integral formula to determine the function 


1 
f(1—#?)$+-c(1— K2#?)!!  (72+-8?)! Pl(s) ds 
¢ t t eo 
Y( ) (1 t2)2(¢? | tT?) } e t 
1 


in terms of the five constants A”, B”, C”, D”, and EB”. These are then 
found from the condition that the singularities in (31) and (32) should 
cancel and from (35). We can then proceed as in section 2 to determine 
the surface pressure distribution. 
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List of symbols 
Throughout the field 
r,y,2 rectangular Cartesian coordinates with origin at the apex of the 
wing, z measured downstream and 2 in the transverse direction (to 
starboard). 
Port jet 
W.p,u values of the undisturbed velocity, density, and Mach angle 


respectively. 


1,v,w components of the additional velocity. 


ty), 4¥;, 7", non-dimensional variables defined by 
v/(2 tan p) Ly r cos 6, 
y/(z tan p) yy rsin @. 


3 angle of sweepback. 
cot u cot B. 
t r/Ssin @+a(1 r=)? cos @}. 


p additional pressure. 


Starboa / jet 


In general, quantities corresponding to those defined in the port jet are 


denoted by the same symbols with a prime added. Exceptions are 


cot uw cot 3. 
t ry’ Kisin @+-i(1—r’*)) cos @!. (N.B. 6 has the same meaning in both 
jets. 


Through Ut the field 
K (tan ; (tan p’). 


( (pW)/(p’W’). 
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\ SPECIAL TYPE OF GROUP DISPLACEMENT FOR 
er mn . , mm mny ’ 
USE IN THE RELAXATION TECHNIQUE 
By R. H. WOOD (Building Research Station, Garston. Herts.) 
[Received 14 December 1950] 
SUMMARY 
The ‘relaxation’ technique, as developed by Southwell, Allen, Fox, and others, is 
now accepted as a powerful method for the solution of many problems in physics and 
engineering. The speed with which the ‘residuals’ can be eliminated depends to a 
large extent on the type of governing equation, and the solution of the biharmonic 
equation in particular can prove to be extremely laborious. In this note a special 
type of group displacement is discussed which can speed up the arithmetical work 
quite considerably. 
1. Special group displacements for the elimination of a single 
residual 
Iv is well known that in the case of the biharmonic equation, unit ‘block’ 


displacements of the type shown in Fig. 1 involve the booking down of a 





Fic. 1. Relaxation pattern for a unit block displacement, 
size 4 4 (biharmonic equation). 
The circled numerals indicate the changes made to 
residuals by the displacements which are shown by the 
numerals underlined. 


large number of residuals, whilst at the same time little is achieved by such 
a process. Whilst searching for group displacements to eliminate certain 
forms of grouped residuals the author was led to the idea of devising special 
group displacements to eliminate a single residual without introducing any 
new residuals at the immediately surrounding points. This process may be 
extended ad lib. according to the size of the field. The method of deriving 
one of these ‘extended patterns’ is given, and the other patterns will then 
appear self-explanatory. 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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Commencing with the smallest of these extended patterns the required 


R displacements x, y, z are written down so as to preserve the necessary 


symmetry (Fig. 2), the aim being to transfer a residual of 100 units at the 
mesh point w to surrounding points without introducing any fresh residuals 


j 
} (0 aN " 
— VY — a 
: Jan) 7 se 
be a nth é 
i y >) 44 r \ 
isle >a 
(a) 0 
Ch 
Ol Fic. 2 isplacements f Fic. 3. Extended pattern for a group- 
ended patter displacement designed to eliminate a 
residual of 100 units. 
| The necessary equations are clearly 
20x0-+-32y— 8z LOO, 
Sx— 25y+- 162 0, 
Pa l6y 22-2 0. 
| for which the solution is 15-506, y = 7-595, « = 4-114 to a sufficient 
order of ace uracy. 

This gives rise to the extended pattern of Fig. 3. In like manner the next 
| two extended patterns are as shown in Figs. 4 and 5. As the field size 
| increases the degree of over-relaxation required becomes quite phenomenal, 
| not only at the centre of the pattern but also at surrounding points. 

Other features of spec lal note are: 

a) All movements are of the same algebraic sign and are symmetrically 
| placed 

sucl b) The largest of the transferred residuals is considerably reduced in 
tall magnitude 

cial | c) Only two rings of residuals have to be booked down as against a 
any | possible five rings or more when using other types of group relaxa- 
\ be tions (e.g. Fig. 1 

ving (d) Under favourable circumstances, by employment of the larger pat- 
thet terns it 1s often possible to transfer part of the initial (100) residual 


over the boundary where residuals do not matter. 
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ing on the conditions. 
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stresses round the boundary provides us (by integration) with values of 
éx/éx and éy/éy around the boundary, where y is the stress function. When 


5. Extended pattern nad aiid of residual of 
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2. Extended patterns near a boundary coinciding with the mesh 


There are several kinds of extended patterns near the boundary depend- 


The extensional problem for a plate subject to known 


29-60-3662 
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Fic. 4. Extended pattern for elimination of 


residual of 100 units—4 x 4 square. 
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the boundary coincides with the mesh points some convenient ‘extended 
patterns’ for points next to the boundary can be derived. These will also 
be applicable to the flexural problem when the slope is given at right angles 


to the boundary. In the extensional problem, having once obtained the 


i 
Woy v4 
~ ~~ ~~" 
4 ao 
an liek 
oe CU 
Fic. 6. Edge pattern near a boundary. 
Fic. 7. Triangular edge pattern near a boundary. 
Fic. 8. Extended pattern—corner case. 


residual for points next to the boundary the movements at such a point are 
wcompanied by similar movements at the fictitious point outside. 

With these conditions it is easy to derive extended patterns of group 
displacements of the rectangular type (Fig. 6) or alternatively of the tri- 
angular type (Fig. 7). These may be enlarged to any size desired. 

In the corner of a rectangular field the movements required are much 
smaller and the first of the extended patterns for the extensional prob- 


lem is shown in Fig. dS. 
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Similar patterns may be evaluated for the flexural problem with a 
straight, simply supported edge by giving the fictitious points next to 
the boundary movements of opposite sign, and reforming the equations 
accordingly. 


3. Recommended procedure in the use of extended group move- 
ments 
The following practical recommendations will be found useful during the 
arithmetical work. 


Fic. 9. The use of the edge pattern (Fig. 6) 
adjacent to a corner. 


|. It is advisable to use the largest possible pattern within the given 
boundaries, so as to transfer as much of the residual as is possible over 
the boundary. 

2. The movements have been quoted to a fair degree of accuracy in order 
to allow large residuals to be eliminated without loss of accuracy. For the 
most part the movements may be evaluated on a slide rule, one setting of 
the slide being usually sufficient for the whole pattern. As the residuals 
decrease in size the displacements do not have to be quoted to the same 
number of significant figures. 

3. Since the rate of convergence is considerably increased (in the author's 
experience some 3—4 times) it is worth while to indicate the outline of the 
pattern of displacements (as in Fig. 9), and if the work is accidentally inter- 
rupted there is no difficulty in tracing the movements which have pre- 
viously been made. 

4. Where a group of positive residuals have to be eliminated the pattern 
as a whole may be ‘over-relaxed’, thus introducing a strong negative resi- 
dual at the centre of the group, leaving the immediate surrounding positive 
residuals unchanged. The whole group is then easily liquidated. 

5. For mesh points next to a centre line of symmetry the extended pat- 
terns may be used taking care to duplicate ‘reflected’ movements as well 
as reflected residuals. Here again the whole pattern is ‘over-relaxed’ to 
counteract any reflected residual at the point of application. 
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6. Whenever the application of an extended pattern causes displacements 

to to be booked down next to the boundary it must be remembered that the 
ns pattern was based on the assumption of a continuous domain, and therefore 
no corresponding displacements of the fictitious points were assumed. In 


that event extra residuals equal and opposite in sign to the displacements 


C- , ‘ 
should be booked at points next to the boundary (in the extensional 
case). The same remarks apply to the use of an ‘edge’ pattern (Fig. 6) when 
applied to a point adjacent to a corner. This operation is shown in Fig. 9. 
4. The analogy with ‘influence’ displacements 
When a very large extended pattern is used to eliminate a single residual 
it is obvious that there will be very little ‘return’ of residual from points 
near the boundary (Figs. 6, 7). Consequently for points in the centre of 
large field, provided the largest possible pattern is adopted, only the 
peripheral displacements undergo any marked change (to suit the boundary 
conditions) whilst the central displacements remain substantially un- 
altered. Evidently the displacements near the centre may be regarded as 
ipproximate influence movements, consequent upon the ‘application’ of a 
single negative residual at the centre of the pattern. In favourable circum- 
stances they are analogous to a deflexion diagram. 
vel r : YY Y/Y i ae 
det 1 
thi “ 
ro Ri * 
1als 
Line wy 2: 
ors T 
tine bU4 (22 
tel 
pre 
| IW-O4 
ter! Fic. 10. Laplace equation—extended pattern for 
resi / squal f side 6 mesh lengths. 


5. Possible uses with other equations 


The procedure, of course, can be adopted for other types of equation, 


pat but such a procedure is obviously unwarranted in the case of the Laplace 
well | equation. However, it often happens that errors are made in the calcula- 
1’ to tion of solitary residuals. There is no faster method of correcting such a 


mistake than by the application of a pattern such as Fig. 10. 
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It is of interest in the table below to compare the largest (central) 
displacements required to transfer a residual of 100 units in the biharmonic 
and plane harmonic cases. 





Btharmonic Plane harmonic 
5 lve - 
Ize Centre l Centre 
(mesh length) | displacement Diff. residual | displacement Diff. 
° 5000 40 25°000 2570 
10°506 12°5 
2 15°506 25°85 37°50 12°5 
14°395 6°73 
4 29°801 222 44°23 77 
18°243 4°006 | 
6 | 48°044 I9°0 48-8096 5°5 Vai 


sys 

The above table shows quite clearly the phenomenal degree of over- = 

. . . 7 . . . n 
relaxation required in the case of the biharmonic equation, and the depen- ae 


dence on the size of the domain. 
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GRADED NETS IN HARMONIC AND BIHARMONIC 


RELAXATION 





By D. N. DE G. ALLEN and 8. C. R. DENNIS 
(Imperial College, London) 


keceived 27 February 1951] 


SUMMARY 


The use of graded relaxation networks was introduced by Allen, Southwell, and 
Vaisey (1) in connexion with the solution of equations of Poisson’s type, but no 
systematic extension of the idea has ever been suggested in applications of the 
relaxation method to solve other partial differential equations in two dimensions. 
In response to several requests the present authors have considered how a graded 
net can be incorporated into the relaxational technique of solution of the biharmonic 


equation ; they came independently to the same conclusions and this paper describes 


how the diffi ulty has been resolved. 


1. Iv is a well-known feature of the relaxation method that its application 
to solve a biharmonic equation is considerably more difficult than its use 
to solve a plane-harmonic or Poisson’s equation; the difficulty is not funda- 
mental but only one of degree in that the convergence of the arithmetical 
process is much slower in the first case than in the second. Consequently 
it would seem to be even more important to develop computational devices 
which will speed up this process for solution of the biharmonic equation 
than it is for solution of a Laplace’s or Poisson’s equation; most of the 
devices in common use—such as over-relaxation, block and group relaxa- 
tion, ete.—were originally propounded in respect of plane-harmonic 
relaxation and subsequently extended to bring biharmonic relaxation 
within their scope. But one most important aid to plane-harmonic relaxa- 
tion has never been extended to apply to biharmonic relaxation; this is the 
use of a graded network, and it is the more remarkable that no efforts have 
been made to extend this device because it would seem to be the one above 
all others which could be put to greatest use in the solution of biharmonic 
problems. 

2. The value of a graded net lies in the fact that it enables the total 
number of nodes in a relaxation network to be kept to a minimum whilst 
allowing a sufficiently fine net to be employed in such parts of a region of 
integration of a given equation as considerations of accuracy may demand. 
A single grading in a net enables us systematically to join up square meshes 
of one size h in one part of the region with meshes of size $h in another 
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part; Fig. 1 shows a part of such a region with a coarse net of side h (to 
the left) connected up to a finer net of side 3h (to the right) by means of 
an intermediate ‘diagonal’ net of side h/\2; this diagram is typical of the 


method of grading as used for plane-harmonic relaxation. The most 
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striking application of its use is in an example given by Southwell and 
Vaisey (2); in their example 2 a net is graded six times from a coarsest 
size of h = 1/2 to a finest size of h 1/128. As a result, their solution of 
the axially symmetric irrotational flow of an incompressible fluid through 
a Borda mouthpiece is able to be presented on a network, of sufficient 
fineness in that comparatively small part of the whole region of flow where 
such fineness is essential for the required accuracy, but which contains in 
all only some 700 nodes. If the finest net had been extended throughout 
the whole region the number of nodes would have been about 1-5 « 10®—a 
number quite impossible to contemplate since for plane-harmonic relaxa- 
tion it may be safely said that it is virtually physically impossible to work 
any solution on a net containing more than about 2,500 nodes. In bi- 
harmonic relaxation, with its noticeably more complicated (or, at least, 
more extensive) relaxation patterns, the maximum number of nodes which 
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to one can contemplate introducing into any solution is considerably less, say 
of of the order of 1,000: thus there is even greater need in biharmonic 
he relaxation for graded nets than there is in plane-harmonic relaxation. 
st 3. Fig. 2 represents a typical group of thirteen nodes of a square net 
_ 
10 
* 
°6 73 5 
bag rs ° 
7 93 %0 1 7 
. > = | 
4 8 
e/2 
Fic. 2. 


of side h. In the solution of a Poisson’s equation 
V2w+ Z(x,y) = 0, (1) 
we define a residual at a typical node 0 by the formula 


F, > w,—4wy +h? Zp, (2) 
where 5 w, stands for the sum of the four w-values at the points 1, 2, 3, 


ind 4. Similarly in the solution of a biharmonic equation 





Viw = Z(x, y), (3) 
and we define a residual at 0 by the formula 
a Fy, = 8 > w,—20w,—2 ¥ w3—)> wythtZp. (4) 
ail It is essential to the relaxation method of solution of any equation that 
‘ent } at each node of the chosen net, where a w-value is to be determined 
sa numerically as part of the solution, a corresponding residual must be 
si te calculable. For Poisson’s equation (1) this is easily seen to be possible at 
out ill nodes in the region of the grading in Fig. 1; for all nodes on the mesh- 
‘ | lines (i), (ii), and (iii) may be regarded as ‘belonging’ to the coarse net and 
-— | so their residuals are defined by (2); all nodes on the lines (vi), (vii), and 
cai | viii) are to be regarded as belonging to the finer net and their residuals 
bi re therefore defined by the formula 
Ast K, > wy 44 | (Sh)?Zo; (5) 


' — 





nodes on the line (v) may be taken in two sets—those like B, E, M, etc., 
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which may be taken to belong to the coarse net with their residuals defined 
by (2), and those like U, V, ete., which may be taken to belong to the 


finer net with their residuals defined by (5). The only remaining nodes are 


those like A on line (iv) and these can be regarded as belonging to an 
intermediate ‘diagonal’ net of side h/\2 inclined at 45° to the sides of the 
coarse and finer nets of sides h and 3h respectively. By virtue of the 
property of the operator V? that it is invariant for any rotational change 
of axes it follows at once that the residual at the node A may be defined 
by the formula 


9 


7 h r, 
F WptWotWpt+uUp teat (5 Z,: (6) 


Hence, for plane-harmonic relaxation, residuals can be defined at all nodes 
of the net in Fig. 1 by means of formulae of one or other of the types (2), 
(5), and (6). 

4. At first sight it might appear that to grade the network in biharmonic 
relaxation would require a much more complicated intermediate region 
than that shown in Fig. 1—for the reason that the standard biharmonic 
residual formula (4) extends over four mesh-lengths, and therefore it might 
be expected that the diagonal net would have to extend at least over four 
mesh-lengths instead of only two, as it does in Fig. 1. However, this is not 
the case, and the difficulty is immediately resolved when it is realized that 
the operator V4 can be written as V?. V? and that these two factors V? can 
be separately replaced by finite-difference operators which are not neces- 
sarily associated with the same size of mesh; the grading shown in Fig. | 
will then be found to be adequate for use in biharmonic relaxation. The 
actual determination of the residual formulae appropriate to the various 
nodes in Fig. 1 will sufficiently explain the method. 

5. In Fig. 1 nodes on the lines (i), (ii), and (iii) may be regarded as 
belonging to the coarse net and their residuals are therefore defined by (4); 
nodes on the lines (vii) and (viii) should be regarded as belonging to the 
finer net and their residuals are then defined by the formula 

F, = 8> w,—20w,)—2 ¥ ws— > wot (4h)*Zp. (7) 
There remains to find residual formulae appropriate to the nodes on lines 
(iv), (v), and (vi). We consider firstly the nodes like A on line (iv); we 
wEne (Viw) , = {V2(V2w)} ,. 
Now firstly replace the outer V* operator by its finite-difference approxima 
tion on the diagonal net: 


V2 


(." ) (v4, (Vw), T (V2w),. T (V2w)p (V2w)_ 4(V*w) ,: 
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ed next replace the remaining V? operators by their finite-difference approxi- 
the mations—at B, C, D, and E on the coarse net and at A on the diagonal 
Are net. We have 

” h?(V?w), Upto tWgtUp—4wp, 

“s h*(Vrw), Wp Wy + Wy +Wp— 4We, 

si h?(V*w)p Ww, WotWetwy 4u D 

ned h?(V*w); Wy+Wpt+Wpt+wy— 4p, 

h 2 
und | : 5) (V2w), Wt wr, Wpyt+Wp—4w,. 


Hence, In all 





‘ m h{—,) (V4w) ' 32w,—10(w_pt+WetWptue_)H 
(Ww, Wy bwWytWytwR ty ty tey), 
mi so that the residual at A appropriate to the biharmonic equation (3) is 
i0l . . 
> F (Wp W, Wh Wr) l6w, 
7 
: h\4, 
ight 9 (Wy Ww, Uy Wy Wye Wy + Wy Wy) | > Z,: (38) 
four ee 
sii Residual formulae at nodes on the lines (v) and (vi) may similarly be 
10 
‘| obtained; on line (v) nodes like B, #, and M can be taken as belonging to 
l 
| : : ; : . 
the coarse net and their residuals are given by the formula (4); at the 
Cal ; . 
alternate nodes on line (v), typified by U, the residual formula is 
ces . . 
ig. | F. 8w, 1 (w, w, bw, Tw, (wptws Wa)— 
The 2(Wetwp)—](wptugtwytwy)+(bhyZy. (9) 
IS ° . ° —— » ° ° 
10u On line (vi) nodes like P and T' belong to the finer net and their residuals 
are defined by the formula (7); at the alternate nodes on line (vi), typified 
d as by Q, the residual formula is 
(4 
ny if o : ” a ‘ Slean | on 
the Fy (3wy+8(wpt+uy, Sw, — 19wg—2(wptwy) 
Wot W, w Wy +wy)—T(wptwpt+uy) (3h)*Zo. (10) 
- Hence biharmonic residuals can be defined at all nodes of the net in Fig. 1 
, by use of formulae of one or other of the types (4), (7), (8), (9), and (10). 
ines : ; : 
From these formulae corresponding relaxation patterns can be straight- 
Wwe . ° 
forwardly constructed in the usual manner. 
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POLYNOMIAL APPROXIMATIONS IN) PLANE ra 
ELASTIC PROBLEMS inte 


the 
' rom Foe . ’ ; s 
By C. A. M. GRAY (The University, Sydney, Australia) ' 
e | 
{Received 17 October 1950] 
SUMMARY If 
, \ 
A general solution of plane elastic problems using the complex variable has been 
developed by Muschelisvili. His solution, for simply connected areas, requires the the 
determination of two anal:-tic functions, given by two integral equations, in which and 
the loaded area in the z-plane is transformed to a unit circle in the ¢-plane. The app 
integrations involved have been carried out only for cases in which the mapping QO’ ; 
function is rational. In this paper a method of solution is derived when the r 
Taylor’s series for the mapping function is known. This leads to solving an infinite 
set of simultaneous equations in which the unknowns are the coefficients in the 
Taylor’s expansion of the analytic functions defined by Muschelisvili. These equa- ; 
tions are readily solved by an iterative process, giving results of any required degree Inte 
of accuracy. To illustrate the application of the method, the problem of a square Q 
with equal and opposite tensions directed along a diagonal, and applied at the 
extremities of that diagonal, is given. 
1. IN reference (1) Muschelisvili (see also 2) has shown that problems in 
two-dimensional elasticity with a given surface loading may be reduced to 
determining two functions Q(z) and w’(z) given by the integral formulae, we 
— SF as : soe OO r ££. ian OF 
QJ f ()] = | (F,+7F) -— — | F (oe) Q, (6) +constant; (1) 
2770 . o—€ 2m, ao—C 
y y 
or S Wiscs cee PoP ae da 
w'| f (f)] _ | (F,—iF,)—; -——. | f (@)Q)(c)—— + constant. (2) 
270 | o—C 27 o—¢ 
y y 
, nena . = whe 
In these formulae z = f(¢) is the relation conformally transforming the 1 
loaded area in the z-plane to a unit circle y: |¢| < 1 in the ¢-plane, o = e” 
is the value of f on y, and bars denote the conjugate functions. Further tion 
p dQ d ees , 
2, = {QF (Oh, (3) 
dz d¢ dz whe 
and F,+-<F, is the transformed value of the load function 
fi+ife = 4i | (X,+i¥,) dp, (t 
. an 
where the line integral in (4) is taken round the loaded area in the z-plane, 
[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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"are the components of the applied surface tractions in the 
| 


and y directions respectively. It is to be noted that the origin of the 


integration in (4) may be chosen quite arbitrarily, since the constant 
thereby introduced can have no effect on the stresses. 
1 it can 


Since z is an analytic function of in the unit circle, y: |¢ 


be expresst d as the powel series 


If we take a finite number of terms, say, in this expansion, and use 
them as an approximation to z, the integrations in (1) can be carried out 
and a value Q' for Q. obtained. But as n is increased, the polynomial 


approximation to z continually improves, and consequently the value of 


Q' will approach the true value 2). Thus, by substituting the polynomial 
ne 
> &¢ 
I 

into (3), Q. can be obtained as closely as desired if n istaken large enough. 


Q bo +6, 6+6,¢ 
aii fF _. do | "9 
(f iF) ” PoTPiStTP2S” 
dC \27 ~ o-—| 





b.4 
Do Ly 109 Uy 01 Us b,, 1U, | Po 
0, U; 2b, u 21D Uy b,, »U, Pi 
(9) 
) 
b_, 2u.b,-» nu,, Do nbo U,, 2 
where the bars denote conjugates as usual, 


/ 


To complete the solution we require w’(z). From (2), with the substitu- 


tionz = ¥ uf and Q, S b. 0, we obtain w’(Z) as 
I 
? a Co ¢,¢ c.¢* 
1 
where - 
, ‘ 7.4 
Co lo 4, 9) iy by+-thg bg | 
Cy q4 1, Dg li, bs lish, | 
Co qo ti, bs lin by li, b; | etc 
| 5 . da oe 
an : (F, I PoT NSS” 
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Thus the solution is completed, and from equation (6) the stresses, x2, ry, 
yy, may be determined 

rutyy = YOQ'(z)+0'(2)] ) 6) 
a ~ —-_~ ” - 2 = ° ) 
ee —YY + vey 20Q"(Z)+0"(2)| J 


In these equations the dashes denote differentiation with regard to z. 
| 




















J 7 
~ * 
| ty 7 
\ | S 
2 y 
7 y 
8 O rd x 
N yy 
f \ 
. % 
\ 
q “i 
. iD 
\/& 2a _ — 
Z PLANE = PLANE 


2. As an example, consider a square, Fig. 1, side 2a loaded at opposite 
corners with a unit tension directed along that diagonal. The transforma- 


tion of a square, side 2a, to the unit circle is given by 


z= (—2ia/k)(C—-105+ -04166729 —-0204¢%...); (7) 
where k 1-85407. Moreover, with this loading, 


Fi+1iF, 4ie’7/4 for ABC, 


0 for C to A. 
l fe : os d 4 im/4 evi7/4__¢ 
Hence : | (F, +7) ae . loz! — =, 
27m , “o—l Qr = Ne8t7/4@_—_ 
Se ieee , ro a 
and a . | (fF, +7) = 4i/a7+40?/7+-4iC4/a—.... (8) 
dl| 277 . “o—C 


Substituting (7) and (8) into (5), and noting, from (8), that: 


(i) the odd b’s are zero, 
(ii) by, bg, etc., are purely imaginary, 
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we obtain the following sets of equations: 


2b) —-1b,+--041667b,—-02404b,,+--0160840D,.... 2k/az, 
by + 1-20834b,—-12020b,+--08042b,.—-058605 ¢... 2k/az, 
6 ~~ ~_ ne] - ~4r 
15b)—+71636b,+- 1-1447b,—-10547b,5+--08122b4.... 2h /a7r, 
62504) + °58410b,—-65235b, + 1-1173b,.—-09382b,.... 2k/arz, 
546866, —-51174b, + -52841b,—-62268b,,+-1-1011b,, = 2k/az, 
-49220b, + -46293b,—-46407b,+ -49985b,. *60527b 4... 2k lax, 
etc., 
und 5 . . — 
1-36, —-125b,.+--07212b,,—-04825b,,+ -03516b,.... 2k/azn, 
‘79167b,+ 1-1683b,—-11259b,,+-08203b,,—-06317),.... —2k/az, 
' -63944b, — -67692b,+ 1-1289b,,—-09426b,,+-07939b,.... = 2k/az, 
‘5537 8b, ‘55079, ‘6 153646 1-1083b,4 ‘O89180,¢... 2k aT, 
-49609), —-48398b,.+--51213b,)—-61296b,,+ 1-09526,,... = 2k/az, 
etc. 
In this case the equations can be solved by iteration. Table 1 gives the 
rs results of five iterations with b’s accepted values and their estimated 
maximum errors. 
TABLE | 
—— 
~ ~ ‘ in 2/ be - b. m7 2h 
413781 347011 
*4309 505 259974 
site 407739 342967 208886 
ma 391824 33134 298357 
S I 54 250051 
76725 "318164 ‘281791 
a ) I (-31 1) 26 I 
Hence Q’ (2k/am)§-465-+--7007122— -3624—-3 1728... (9) 
ind w’(z) = ¢y+(4€/7)|[ 1-7492+-73761f? — -550004 — -4405028... |. (10) 
A graphical determination of these stresses has been obtained by 





} Frocht (3). In Table 2 values of the principal stresses along the diagonals 


| 
| 
| 
| 
j 


re given as computed for (9) and (10), and as determined by Frocht. 


TABLE : 
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From the table there is a discrepancy of about 5 per cent. between the 
two sets of values at the origin. However, Frocht’s values are too high 
as he gives 1-05 for the sum of the normal stresses along OU’; the values 
from equations (9) and (10) being 1-02. Further, if we use the errors 
estimated from the iteration process, limits can be set to the values 
obtained from (9) and (10), which are, for this case, approximately 

-012/a. This gives an extreme error of 1-5 per cent. in the values of the 
major stress at the origin. 

Generally the method is available wherever a power series expansion is 
known for the transformation z = f(¢). In the case of polygonal boun- 
daries such developments can be constructed from the Schwartz-Christoffel 
theorem, although the convergence of the iteration process may not occur 
without modification of the equations. However, preliminary solutions 
for a long rectangle with normal loads on the long sides have shown that, 
in this case also, the convergence is sufficiently rapid to be of practical 
importance. 
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DIFFERENTIATION OF FOURIER SERIES IN STRESS 
SOLUTIONS FOR RECTANGULAR PLATES 
By A. M. WINSLOW (University of Washington, Seattle) 
Received 16 January 1951] 


SUMMARY 


Certain mathematical conditions of differentiation are necessarily included in a 
Fourier solution which is based on the hypothesis that a stress function and its 
partial derivatives can be represented by Fourier series. To illustrate these conditions, 


a summary of essential theory is followed by examples, using full-range Fourier series 


und also half-range Fourier se1 


In these examples of stress solutions it is found that the assigned boundary trac- 
tions provide equations sufficient to determine the Fourier coefficients. Additional 
mathematical conditions, however, are imposed by the hypothesis of termwise 
differentiation. Unless it can be proved that these additional conditions are fulfilled, 
the conclusion is that the hypothetical representation by Fourier series may not have 
sufficient generality to satisfy all required conditions and furnish a solution. 
Notations 


THE following notations are frequently used: 


sa maa. 
, nab. 
2a = length of plate in x-direction, in full-range example. 
2b length of plate in y-direction, in full-range example. 


a,b = lengths of plate in (2, y)-directions, respectively, in half-range 
example. 


Other notations are conveniently defined where introduced. 


Differentiation of full-range Fourier series 

Essential principles of differentiation, expounded by Hobson (1), are as 
follows. We consider an even function F(x) and its derivative F’(x), which 
comply with sufficient conditions to permit representation by Fourier series 
in the interval (—a,a), 


x 


0 . — a 

F (a) ~+ > a,,cosB,, 2, (1) 
a mt 1 

F(x) = > 6b, sin B,, x. (2) 
m=1 


If F(x) is a bounded function and is piecewise continuous in (—a,a), 


- tN ; . 
b BnQm-+— > {F(xg—0)—F(xg+0)}sin B,, ta, (3) 
ee 
the summation § referring to the points x, of ordinary discontinuity of 
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F(a) in the interior of (—a,a). This value of b', in (3) is obtained by 
applying integration by parts to the expression 
a 
if 
a, = F(x)cos B,, x da. (4) 
—a 


In a similar manner, with an odd function f(x) 


x 


f(x) > Dn sin 8B, x, (5) 
m=1 
ial %, NS» 
J (*)= s+ > Am COS B,, ®, (6) 
2 m=1 


it is shown that 


a = 1 Fa 0)—f(—a+0) +251 f(ta—9) —f (ta +9)}, (7) 


+-[(—1)™{f(a—9)—f(—a+0)}+ & {f(ea—) —f (wa + O)jeos B,, xa]. (8) 


‘The expression (7) is found by oe 
a = 2 f'(x) da. (9) 


If there are no points of iia xq in the interior of (—a, a), 


’ 


ay = ‘[fa—0) f(—a+0)], (10) 


An Bn bm + —" [f(a —0)—f(—a- 0)]. (11) 


a 


Hence, if f(a) is an odd function represented by the Fourier series (5), 
termwise ison of the Fourier series does not represent f’(x) unless 


f(a—0) = f(—a+0) = 0. (12) 

These results are exemplified by f(x) = 2°, f’(2) = 327. Throughout 
(—a,a), termwise differentiation of the Fourier series for x? does not 
represent the derivative 32°; the derivative is expressed by (6), (10), 
(od). 

Relative to a function f(x) which is either odd or even, in order that 
termwise differentiation of the Fourier series for f(x) may represent the 
derivative f'(x), the necessary condition defined by (3), (7), (8) can be 
stated in general terms. The periodic function represented by the Fourier 


series for f(z) must be continuous for all values of 2, both within (—a,4) 
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and at the ends of the interval, x ta. In particular, if f(x) is an odd 
function, f(x) must converge to the value zero at x +a; that is, 


f(a 0) f(—a t()) 0. 
This condition of complete continuity is inherent in every hypothesis of 
termwise differentiation, and must be included in the mathematical 


solution. 


Example using full-range Fourier series 


With coplanar stresses, using Airy’s stress function ¢(z, y), 


o*h ond od ‘ 
Oy >? Oy o..9? Try . ‘ : ? (13) 
oy? Ca* cucy 
where ¢ also satisfies the compatibility equation 
ad a4 4 
C C C 
V‘d = 0, Vi = —+2—___4 ; (14) 
cat ox*oy* © cy? 


When ¢(x, y) is an x-even and y-even function, a typical expression for 


d 1S Y y = 
ri) P S,+4g, (15) 
where P, a polynomial satisfying (14), is 
1 4 2972 2 2 
.¥ , 2 ; —), 2 ee ' 
P — C,2+C, (Cy +C,) 7 +OF +O,—, (16) 
24 “24 7 2 2 
S, > S(Y)mCcosB,, x, (17) 
m l 
S, = ¥ g(x), cosa, y, (18) 
n 1 
f(Y)m = F.,BmysinhB,, y+ F,, cosh B,, y, (19) 
g(x), Gx, «sinh a, 2+ G, cosh «, x. (20) 


In the series S, and S,, each separate term satisfies the compatibility 
equation (14). Therefore an inherent hypothesis is that successive deriva- 
tives of ¢ used in (14) are correctly represented by termwise differentiation 
of (15). Also, it is inherently assumed that, for successive derivatives of ¢ 
up to 644/éx2éy?, the relative order of differentiation is interchangeable. 


As an example, consider a rectangular plate having boundaries x La, 
y -b, with assigned coplanar tractions, 
CO, A+ By?, T. 0, at x ra | 
‘ ry } ; (21) 
CO, Q, Ti 0, at y + bS 


A method of solving this problem has been described by Lardy (2) and 
Pickett (3) in terms of (15), (16), (17), (18), using P = Cy?. More generality 


, 


Cy2+C'x. (16’) 


would be prov ided by P 
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With this slight revision, the solution will be outlined here, in order to 
explain additional conditions of termwise differentiation. 

Referring to (19), (20), two of the boundary conditions (21) determine 
the coefficients F’,, in terms of F,,, and G; 


» in terms of G,. Then the other 
two boundary conditions (21) give a solution for F,, and G,,. First, con- 
sidering shearing stresses, and using (13), (15), (16’), 


x 


9 78 
o*” ‘ y ° , . 
‘ey = “3 ie 3 J (Ym Bn sin B,, aa Z g (2 )n %, SIN a, Y- (22) 
OX coy m=1 n=1 
In (22), as a necessary condition for further termwise differentiation with 


respect to x, the Fourier series in sin f,,« must represent a function which 


m 


converges to zero at x +a. Hence the assigned condition (21), 7,,, = 0 
at x +a, is satisfied if g(a) = 0, giving 

G,, = —(1+.a, acoth a, a)G,. (23) 
Similarly the condition 7,, = 0 at y = +5 is satisfied if 

Fn -(1+8,, bcothB,, b)F,. (24) 


Results (23), (24) permit representation of S, and S, as linear expressions 
in F,, and G@,. Applying the condition (21), ¢, = A+ By? at x +a, and 
using (13), (15), (16), 


5 
A+ By? = 20+ ¥ (—1)"f"(y)n— > g(a), 2 cos x, y. 25) 
m=1 n=1 

An expression for C is obtained by the usual Fourier method of multiplying 
each term by dy and integrating from y = —b to y= b. Im (25) the 

Fourier coefficient of cos «,, y is 

b 
g(a), 2% = ; | |- By?+ > ( 1)"f"(Y)m|eos x, y dy. (26) 
4 n 


1=1 
b 


Similarly the boundary condition, c, = 0 at y = +5, gives an expression 
for C’, and also 


life ‘ i 
F(9)m Bin . | > (—1)"g"(x),, cos B,, x dx. (27) 
“_ t=1 


Each term of the integrals in (26), (27) can be evaluated by two succes- 
sive integrations by parts. It is seen that (26), (27) are homogeneous linear 
equations in B, F,,, G,,, which can be reduced to the form 


m? n 


n m,n m? 


G, = KE, B+ > CF 
1 


m 


F b Ae (27') 


m 


nym n? 


n=1 
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where C.,,,,, Cy A, are known constants. If finite summations up to r 
terms are used, there are 27 simultaneous equations (26’), (27’) in the 2r 
coefficients F,, G,. In this manner an approximate solution may be 
obtained, which may converge towards an exact solution when the number 
of terms r is increased. 

This is the solution described by Lardy and Pickett. However, not all 
necessary conditions of termwise differentiation have been applied. An 
additional condition is the hypothesis that termwise differentiation of the 
Fourier expression for ¢°¢/cy* represents the derivative ¢4¢/éy*. Equation 
(25) represents the value of co, = 6*4/éy” at the boundary x =a. The 
evaluation of (26) to obtain (26’) involves the value of 0¢/éy*? on the 
boundary 2 a, and at the ends of the interval y = +6. The following 
3 


expression for 6*4/éy? on x = a is found by operating with @/éy on (25). 


Rearranging terms, 
> 9(@2), G sina, y = 2By — ¥ (—1)"f"(Y) m- (28) 
n l m=1 
In (28) a necessary condition for further termwise differentiation with 
respect to y requires that the series in sina, y must represent a function 


which converges to zero at y + b, giving 
2Bb— > (—1)"f'""(b),, = 9. (29) 
fit 1 


In a similar manner, from the value of ¢°¢/éx3 on the boundary y = }, it 
is found that 
> (—1)"9"(a),, = 9. (30) 


l 
These two necessary conditions (29), (30) must be satisfied, in addition 
to the 27 simultaneous equations (26’), (27’). The manner in which the 
conditions (29), (30) affect the stress solution is shown by applying two 


successive integrations by parts to (26), (27), giving 


' ? ] n a om 
Q\a),, xy > 2 Bb > ( 1)”f (b)n| - 
b Xn m=1 : 
b x 
l as “ on 
hy S ( l nf " (Y)m COs Xn y dy, (31 ) 
iy J oe 
, m 1 
‘ ; ?( } = “ 
f(b), B : l)"g"(a),,— 


1)"g'*(x), cosB,,2 dx. (32) 
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Also, it is found that, if conditions (29), (30) are omitted, equations (26’), 
(27’) can be derived identically either from (31), (32) or from the previous 
equations (26), (27). 

The purpose of the present example is to point out the additional 
requirements imposed on the solution by the conditions (29), (30). Later, 
a similar solution in terms of half-range Fourier series will be derived more 
completely. 


Differentiation of half-range Fourier series 

In some solutions a function is represented by Fourier series in an interval 
of length a which is one-half of the full interval of periodicity 2a. A typical 
representation is shown in Fig. 1, where f(a) and its derivatives are alter- 


{ 
_ 
-—-— 

















an . 0 
= a ae ae 
(a) ro. 1. (b) 


nately even and odd functions in (—a/2,a/2) with the origin at O. This 
a-even function f(x), continuous in (—a/2,a/2), is shown as a half-range 
representation of an x’-odd function F(x’) in the full interval (—a,a) with 
the origin at O’, 


F(x’) > by, sin B,,, 2” (33) 
m=1,3,5 
2, 
db}, . F(x’ )sin B., a’ dx’. (34) 
a 
0 


The derivative F’(x’) is not represented by termwise differentiation of 
the series (33) for F(x’), because F(x’) has an ordinary discontinuity at 
x 0, and also because F(a—0) 4 F(—a-+-0) 4 0. As shown previously 
in (5), (7), (8), in order to permit termwise differentiation, the periodic 
eantion represented by the Fourier series (33) must have complete con- 


tinuity, at 2’ = 0 and 2’ La. 
Changing the origin from 0’ to O, by the relation 2’ = x+<a/2, and 
putting 6b”, = a,, sin(m7z/2), the expressions (33), (34) become 
J (x) Zz Am cos B..2, (35) 
m=1,3,5... 
a/2 


f(x)cos B,, x dx, (36) 


m 
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ay which is the usual half-range representation of the even function f(x). 
us The coefficients a,, are zero for even values of m, including m = 0. 


For the derivative 


al f'(x) > 4b, sinB,, x, (37) 
T, m = 1,3,5... 
a/j2 
re ‘ / | | 
bh. | f ‘(x)sin Bn x dx. (38) 
a 


conditions of differentiation could be found from (5), (8), (33), but are 











al obtained more directly by applying integration by parts to (36), 
al - — 
7 2 [ 2 862 6f , 
f(x)cos B,, x dx | — f(x)sin B,,, x| — f'(x)sin B,, « dx. 
a. [a B }—a/2 ap, . 
a/2 a/2 
> - 
Hence, i _ f(s 0) 1A(—$4 y sin ar Se. 
Pm - - = Bm 
Since 
(5 0] i : 0), Bin = —On Bn +a1(5—9)sin — wm 
| Thus, termwise differentiation of the Fourier series (35) for f(x) does not 
his represent the derivative f’(2) unless 
Be f{q 0): f(- ool 0) 0, (40) 
ith 2 . 
- Referring to Fig. 1 (a), it is seen that this necessary condition (40) for 
ga termwise differentiation requires complete continuity of the periodic func- 
tion represented by the Fourier series (35). 
34) With an odd function /f’(x#), continuous in (—a/2,a/2), and represented 
by (37), a derivation similar to that for (39) shows that termwise differentia- 
1 of tion of the Fourier series (37) represents the derivative f"(~). Referring to 
at Fig. 1(b), it is seen that the necessary condition of complete continuity of 
sly the periodic representation (37) is fulfilled. 
die | These results are exemplified by f(x) = a2. Termwise differentiation 
on of the Fourier series for 2* in (—a/2,a/2) does not represent f’(x) 473: 
the expression (39) must be used. However, with f’(x) 4x3 expressed 
and | by (37), (38), termwise differentiation of (37) correctly represents 
f"(w) = 122%. 
(35 | Example using half-range Fourier series 
The previous example of a rectangular plate encountered difficulty in 
(36 the stress solution, caused by conditions (29), (30) of termwise differentia- 


tion. As a preferable solution of this problem, the following method has 
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been suggested to the author. While this proposed solution is found to be 
not free from difficulties of differentiation, it illustrates applications of 
half-range Fourier series. The assigned boundary tractions are 


o, = A+ By?, Try = 0, at z +a/2 | 
= (41) 
Cy 0, Try 0, at Yy so b/2 J 
With the stress function ¢,(2, y) satisfying V4d, = 0, let 
$1 = bo—$, (42) 
Ay?  By' 
where dy = +: (43) 


The stress components (13), using (43), satisfy all boundary tractions 
41), and 
aaa Vibo = 2B. (44) 
Therefore, a solution of the problem is obtained if ¢ can be determined 
to satisfy zero boundary tractions, and also 
V4d = 2B = p/D. (45) 
By putting the constant 2B = p/D this determination of ¢ is recognized 
as that appropriate to the problem of a clamped rectangular plate, where 
the constant p is the intensity of uniform transverse loading, and the 
constant D is the flexural rigidity of the plate. A typical solution for the 
clamped plate in terms of half-range Fourier series is given by Wojtaszak 
(4). It is essential to outline this solution briefly, in order to show con- 
ditions of termwise differentiation. In accordance with (15), 
¢ = P+8,+8,, 
where the series S, and S,, using odd values of m and n, are expressed by 
(17), (18), (19), (20), and 


> Pp a? . ae P 46 
site nar ' \(; r’). (46) 


First, considering the boundary traction o,, 


op p ( ; 
Cr = = — Ye} 


oy” ~ 4D\4 
Y Sf" (Y)mCOSBy x Y g(x), GB cosa, y. (47) 
m=1,3,5... n= 1,3,5... 


In (47) the Fourier series in cos 8,, x must represent a function which con- 


m 


verges to zero at x ta/2, since this is a necessary condition (40) for 
further termwise differentiation with respect to x. The boundary con- 
dition, oc, = 0 at x = +a/2, is satisfied by (47) if g(a/2), = 0. Hence, 
using (20), 


' — G,(a 


n 


, a/2)tanh(«,, a/2). (48) 
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In a similar manner, to satisfy the boundary condition, ¢, = Oat y = +b/2, 
Fh FAB, 6/2)tanh(,,, 6/2). (49) 
Next, considering shearing tractions on the boundaries, 
c “h pe uU 
Cx y 2D 
o 
> J \Y)m 6, sin B,, x > g (x), Xn sin Xn y. (50) 
m= 1,3,5... n=1,3,5... 
In (50), applying the boundary condition 7,, = 0 at x = +a/2, 
pay ‘ > ° S , : 5 
| S 1 (Y)m Pm Sin(m7/2)-4 > g (4/2), «,s8iMa, y 0. (51) 
i? wnt * n=1,3,5.... nit Ps 


In (51) the Fourier coefficient of sin x, y is 


b/2 
= pay J - c ‘ ~ 
g (4/2), 0 Sf’ (Ym Bm Sin(m7/2)} sin x, y dy. (52) 
b 4D am-fas.. ~ ” si 


fs 


Each term of the integrals in (52) can be evaluated by two successive 
integrations by parts. Also,.the final equation is simplified if F,, and G, 


m 


are replaced by dimensionless coefficients A,, and B,, defined by the 





m 
expressions 
A,, pa (53) 
8 DB? cosh(,,, 6 y 
B, pb " 
G = . (54) 
8Da8 cosh(«, a/2) 
Finally, by using (48), (49), equation (52) is reduced to the form 
’ f o2 \2. sal 
N. B.- S Am| a ” —) sin(mz/2) 3 (55) 
a Din T &y, 
l ) 
. nab |nza/b-+-sinh(nza/b)] . a 
where N,, ) sin(nz/2). (56) 
16a | cosh?(nza/2b) 
Similarly, in (50), applying the boundary condition r,,, = Oaty Lb/2, 
, pe 2 
MA SY B, | — ) sin(nz/2) = 1, (57) 
735 Bin T Xn 
mralmrb/at-sinh(mzb/a)] . a: 
where M,, ~ sin(mz7/2). (58) 
165 cosh?(mzb/2a) 


Equations (55), (57) are Wojtaszak’s final results. With finite summa- 
tions up to m rn r, there are r+-1 simultaneous equations (55), (57) 
in the r+-1 coefficients A,,, B,. However, there are conditions of termwise 
differentiation which have not yet been applied in the derivation of (55) 
from (52). This derivation, using integration by parts, involves the value 
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of 24/éxéy? on the boundary x = a/2, and at the ends of the interval 
y -b/2, that is, at the point x = a/2, y = b/2. 
Using the general expression (15), 
ad @eP aS, Aas, _ 
roe teal oe ee ee ee (59) 
Oxdy* dxdy* dxoy* dcxdoy* 
At the point x = a/2, y = b/2, 
a3, nae 
as - = — : cy = 0. 
oxoy* oy 
Using (18), €2S,/ey? is a Fourier series in cosa, y representing a function 
Lo] e « Ne =) 
which converges to zero on the boundary y = b/2. This is a necessary 
condition for further termwise differentiation with respect to y. Therefore, 
the x-derivative of this function is zero on the boundary y = 6/2 and, in 
particular, at the point x = a/2, y = b/2. Thus, at this point, (59) becomes 


PF —  ¥_—f"(b/2)mBp 8in(mn/2) = 0, (60) 
1,< 


4D m 13,5... 


which, by using (49), (53), reduces to 


A,, sin(m7/2) = 0. (61) 


Ms 


_ 
i) 


3355... 


m 


Similarly, considering the derivation of (57), 


Ms 


B,, sin(nz/2) = 0. (62) 


Ul 


w 


n=1,3,5... 

The necessary conditions (61), (62) must be satisfied, together with the 
simultaneous equations (55), (57). The manner in which the stress solution 
is affected by condition (60) is seen by applying two successive integrations 
by parts to (52), giving 





P 4 [pa x oe - . 
g'(a/2),, x, - Sf" (b/2) By Sin(m7/2) | sin(m7/2)-4 
< t h x2 4D —_ = e m m 
n m= 1,3,5... 
b/2 
9 Pf < 
, s ™ ee 66 
+75 Sf” (Ym Bm Sin(m7/2)sin x, y dy. (63) 
we big M=1,3,5... 


Furthermore, it is found that, if condition (60) is omitted, equation (55) can 
be derived identically either from (63) or from the previous equation (52). 
For a square plate, where (55) and (57) are identical, and 


b= a, B re x g.., 


n n n 
Wojtaszak used a finite summation up to m = r, with r = 27, and com- 
puted values of the coefficients A,,. Also, the author has made finite 


solutions, using 7 15 andr = 47, in order to compare values of A,,, with 
reference to the necessary condition (61). The results are shown in Table 1. 
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TABLE | 


Values of A,,, using finite series up tom = r, for various values of r 
’ 47 r 27 r 15 
4135 1°4135 1°4135 
143 0°0943 0°0943 
90 o*1096 O° 1097 
765 00705 0°0707 
523 00523 0°0520 
05 0°0305 0°0370 
0°0260 0°0206 
158 00188 2°O197 
135 0°01 35 
I 102 o"o1o02 
076 0°0076 
5 0°0057 
O41 0°0042 
>30 0°0031 
1 
14 
»Q 
I 
I 
With r = 15, 27, 47, the successive values of the left side of equation (61) 


are 0-0028, 0-0448, 0-0472. While not establishing proof relative to infinite 
series, these results show no convergence towards agreement with condition 


(61), when the finite number of terms is increased, fromr = l5uptor = 47. 


Conclusion 

The two preceding methods of solution, in full-range and half-range 
Fourier series, have been used to illustrate conditions of differentiation, 
and are typical examples of the following situation. 

When a solution is based on the hypothesis that Airy’s stress function 
d(x, y) and its derivatives can be represented by Fourier series, the succes- 
sive partial derivatives being obtained by termwise differentiation, it is 
possible that the hypothetical representation does not have sufficient 
generality to satisfy all required conditions. The generality of a function 
represented by Fourier series is greatly restricted by the hypothesis that 
termwise differentiation represents the derivative of the function. There- 
fore, in a solution which must comply with assigned boundary tractions, 
requirements of termwise differentiation must be considered as additional 


independent restrictions, unless proved otherwise. 
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In a specific example applying termwise differentiation, it may be diffi- 
cult to prove that the mathematical requirements are satisfied, or that they 
are not satisfied. Nevertheless, whatever mathematical difficulties are in- 
volved, all necessary conditions of differentiation must be included in the 


solution. Otherwise. there is no assurance that successive derivatives of 


the stress function ¢ are correctly represented by termwise differentiation, 
and there is no guarantee that ¢ satisfies the compatibility equation 


Vd = 0. 
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HEAT CONDUCTION IN THE REGION BOUNDED 

INTERNALLY BY AN ELLIPTICAL CYLINDER AND 

AN ANALOGOUS PROBLEM IN ATMOSPHERIC 
DIFFUSION 


By C. J. TRANTER (Military College of Science, Shrivenham) 
Received 20 March 1951] 


SUMMARY 


\ formal solution is given for the flow of heat in the region bounded internally by 
an elliptical cylinder. A particular case of the solution can be used in the problem 
of atmospheric diffusion over a saturated surface in the form of a long rectangular 
strip when wind speed and diffusivity are constant with height. 


1. Introduction 
THE flow of heat inside a long elliptical cylinder with zero initial and 
constant surface temperature has been discussed by McLachlan (1). A 
solution to the corresponding problem for the region outside an elliptical 
cylinder does not seem to have been given so far and a formal solution 
is obtained in the present note. For the internal problem, McLachlan was 
able to obtain the solution in the form of a double series of products of 
the Mathieu functions Ce, ce and to determine the coefficients by an 
orthogonality theorem. There seems to be no corresponding orthogonality 
theorem for the external problem and the analysis employed here is rather 
different 

The solution of the heat conduction problem can be used in the analogous 
problem of atmospheric diffusion over a saturated surface in the form of 
a long rectangular strip whose long edges are parallel to the wind direction 
when the mean wind velocity and diffusivities in vertical and cross-wind 


directions are assumed constant with height. 


2. Heat conduction in the region bounded internally by an elliptical 
cylinder 
If we take the inter-focal line of the cylinder to be of length 2h and use 


the usual elliptical coordinates (€, 7) [(1), p. 170], the equation for solution 


eV ei h2 oV . 
=| (cosh 2&— cos 2) (E Pe 0,0 < 7» < 2n7,t > 9) 
0&< 0n* 2K ot 


where V is the temperature at time ¢, « the diffusivity, and &, is the value 


(Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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of € on the bounding ellipse. Assuming zero initial temperature and 
constant temperature V, on the surface € = &), equation (1) has to be 


solved with 


V>0, &<&E<wO, 0O<yn< 22, ast—> +9, (2) 
V>Kh, 0<n< 22, t>0, E> &, (3) 
V>0, O<yn< 227, t>0, E>. (4) 


Write V as the Laplace transform of V, i.e. 
V—[Vedt (p>), (5) 
0 


then integration by parts and use of (2) gives 


| ot e—Pt dt pV. 
| at 

0 
Multiplying (1), (3), and (4) by e~” and integrating with respect to ¢ from 
0 to oO gives 


ev @ 


— 2q(cosh 2 — cos 2n)V (&<€<0, O<7< 27), (6) 

0€* — @n* 
where q = ph?/4k, (7) 
and V>Vip (0< 7 < 2z, as > &), (8) 
V>oO (0 < » < 2z, &—> 00). (9) 


In view of (9) and the symmetry about both axes of the ellipse, the 
appropriate solution of (6) is 


V 2 ( Qn Fek,,,(&, J)CCg,,( > q); ( 10) 
n=0 
where the Mathicu functions ce,,,, Fek,,, are those defined by McLachlan 


on pages 21, 165 of his book.* The constants C,,, have now to be deter- 


mined from the remaining boundary condition (8), i.e. 


K/p = > C,, Feks,(&, —q)ceen(n, —g) (0 < 9 < 2z). (11) 
n=0 
Multiplying (11) by ce,,,(y, —q), integrating with respect to 7 between 0 


and 27, using (1) pp. 23, 24, viz. 


27 
| ces,,(7, —G)Ceam(y, —g) dn = 0 (m #7), 
0 
27 
9 
ces,,(, —q) dy 7, 
0 
* McLachlan (p. 257) indicates that the Bessel function product series given on pp. 246 


and 248 have some advantages over the series given on p. 165. 
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. | and (1) p. 21, which gives 

be 2a 0 

| cey,,(n, —q) dy (—1)" | > (— 1)" AZ” cos 2rn dn (—1)"27AQ™, 
) 


we have 

















3) ' 2, ony Heke, (€, —@) 

J > 1)? A ___* ce,,,(n, —q). (12) 
4) pP rea’ ‘ Fek,,, (£5, q) r " ) 

The Mellin inversion theorem now gives 
(5) V, t ; 4 ony FeKe,, (€, —p) dx : 
} : (—1) eX AGh __ ce,,,(7, —u) —, 13 
71 yaar J . Fek,,, (& , LL) 2n( I v) A \ ) 
|B 
ym 
(6) A 
(4) 
ss F E 
(9) ~T\ 
the Cc ew, 
10) 
lan 
el 
11) 
nv 
A 





where » = Ah®/4« and the argument of A@” is uy. There is a branch point 


at the origin and, using the method of Carslaw and Jaeger (2a), the line 
integral is replaced by a real infinite integral, derived from the integrals 
- along CD and EF together with a contribution from the small circle about 
— the origin (Fig. 1). 


| 
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The integral round the small circle gives 77 in the limit as its radius 
tends to zero. On EF we put A ve'” and since 


K,,.(ze§**) Lrie7'| J,.(z)—iY5,(z) |, 
use_of equations (9), (1) p. 165, (1) and (15), (1) p. 159, gives 
Fek,,,(€, —p) bi Ce, (€, xh? /4x) + Fey,,, (€, ah?/4x) |, 


the Mathieu functions Ce,,,, Fey,,, being defined on p. 159 of McLachlan’s 
book. Hence the contribution from FF is, writing ah?/4« = v, 





sy ( l * P Seve n2 | yn 4(2n) a §. e anit; v)-4 F e ad walt, Vv ) ce (7 yy” 

. —~ 2n\")> . 

a, m , i Ces, (Ey, v) + Feys,, (Eo: v) v 
0 


The contribution from CD gives minus the conjugate of this. Combining 
these results, we have finally 


V 1.2 s [ ¢ Axvt/h® 4(2n) < 
My 7 
n=0 9 
Ce,, (€,) Fey, (£9. ¥) —Feys, (, v) Cea, (Eo. » by 
San(6,¥) £¥ anl€or) oYan\s p)LeanlSor¥) ce,,,(1,v) ‘ “, (14) 
Ce3,, (Eo, v) + Fey3,, (So, ») . v 


the A\?” being functions of v. 


3. Transition to circular cylinder 

When the ellipse degenerates into a circle of radius a and the confocal 
ellipses become concentric circles, a typical radius being r, h—> 0, € and 
&,—> © in such a way that shef > r, the&—>a. Under these conditions, 
(1) pp. 367, 368, 


A@) +0. (n ~ 0), Am > 2-4, C€o(n, ah?/4x) > 2 


Ce ate vh?/4«) > Don Jun( I(r) Fey,,,( (é, xh? 4K) > Pan Yon(/(2)) 
K K 


where p,,, is independent of 7, and with similar expressions for the functions 
of argument &, except that a now replaces r. Making these substitutions 
and writing a = xu?, we find 


V +i. 2 : cot moter —Y,(ur)Jy(wa)|] du 
% . J?(wa) -Y2(ua) u?’ 


a result in agreement with that given by Carslaw and Jaeger (25). 
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HEAT CONDUCTION 


4. Atmospheric diffusion from a long rectangular strip 

We take the saturated strip (width 2h) to lie in the plane z = 0, the 
longer sides being parallel to the wind direction which is taken to lie along 
the v-axis. The z-axis is taken vertically upwards and the y-axis at right 
angles to the wind. If U is the mean wind velocity, K,, K, the coefficients 
of diffusivity in the y and z directions, and if we neglect the coefficient of 
diffusivity parallel to the wind, the equation satisfied by the concentration 
(x,y,z) of vapour is, under steady state conditions (3), 


eV of, eV of, eV . 
l kK, +.—| K, —}, (15) 
CX Cy, ~ cy CZ| ~ 
with boundary conditions 
] 0, 2z 0, asx>-+0 
V+ r>o0 h<y<h, z>0 ' 
0 . e : (16) 
lim; K, ~0O, x >, h>y>h 
0 | ” O2) 
and V tending to zero at large distances from the strip. 
Write 
x Ut, y K>} h cosh € cos », z = Kihsinhésin n, (17) 
then (15) and (16) become, if U’, = K, are assumed constant, 
Vy a h? V 
tated ~ (cosh 2—cos2y)— (0<E<a0,0<y <2), (18) 
0g* on- yA ot 
with | 0 0 é x 0 <7 t~> 0 
| > Vi, U 1) 7 t > Q, € > 0, : (19) 
. {av 
lim | a, €>8; One < oe 
n—-0,7\ C7 


The last condition is equivalent to assuming symmetry about the major 
axes of the confocal ellipses of which the saturated strip is the limiting case. 
Hence the solution to this problem is identical with that of para. 2 when 


we write €&, = 0, x | and interpret €, , ¢ in accordance with (17).* 

* The result given by Davies is quite different from that given here. It seems that the 
orthogonality theorem used by him in obtaining the coefficients in his equation (3.12) is not 
alid. This defect appears to arise from the assumed vanishing of the first integral on the 
left of his equ ition (4.25). Although as B-->o, Fey.,(B) » 0, its derivative with respect to B 
8 infinite and the product of the function and its derivative is finite as B + 00. 
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SUMMARY 

The aim of this paper is to obtain an estimate of the effect of finite width on the 
rate of evaporation from a plane, rectangular, saturated area into a fully turbulent 
air stream. A rigorous solution for this problem has been previously shown to 
involve Mathieu functions of a type which have not yet been computed; in order to 
obtain an estimate of the effect of finite width a simplified version of the problem is 
discussed. 

The problem considered is that of turbulent diffusion of matter from a plane area 
of infinite length and finite width, composed of infinite line sources, when the mean 
wind direction is parallel to the sources. The strengths of the sources are chosen to 
be constant along the mean wind direction and variable across wind, in such a way 
that the concentration of vapour is constant at all points on the area. The use of 
this idealized model is justified by comparing the theoretical law obtained, for 
dependence of rate of evaporation on width, with experimental results. 

The analysis leads to the theoretical result that the rate of evaporation per unit 
length of the area is not directly proportional to the width, as assumed by previous 
investigators, but to (width)® where, for example, in the particular case of adiabatic 
conditions a = 0-86; the theory also indicates that the index « is independent of the 
mean wind speed and depends only on the degree of turbulence present in the air 
stream. Experimental results are given which suggest that, in wind-tunnel condi- 
tions, a law of this kind expressing the dependence of rate of evaporation on width 


is true at all points along an evaporating area, except very near the leading edge. 


1. Introduction 
ALTHOUGH the problem of turbulent dispersion of vapour above a finite 
evaporating area has been worked out in detail by D. R. Davies (1), when 
the area is bounded by a parabolic curve, it has not been possible in the 
past to derive formulae which would lead to a quantitative assessment of 
the effect of finite width on the diffusion process. It has been previously 
shown (2) that a rigorous solution for the problem of diffusion from a 
rectangular area involves Mathieu functions of a type which have not yet 
been computed, and in order to obtain an estimate of the effect of finite 
width on the rate of evaporation a simplified version of the problem is 
discussed in this paper. 

The problem considered is that of turbulent diffusion of vapour from 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 4 (1951)] 
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a plane, rectangular area, of infinite length and finite width, composed of 
infinite line sources, when the mean wind direction is parallel to the sources. 
The strengths of the sources are chosen to be constant along the mean wind 
direction and variable across wind, in such a way that the concentration 
of vapour is constant on the area. The use of this idealized model for 
making a theoretical prediction of the influence of finite width on rate 
of evaporation is justified by comparing the theoretical law of dependence 
on width with experimental results obtained in a wind tunnel. 

This problem is solved, firstly by applying the source method, used by 
K. L. Calder (3) in another context, and secondly by obtaining a suitable 
solution of the diffusion equation. The approach to the problem following 
the method of sources enables an empirical formula to be obtained for the 
factor a, in the power law of eddy diffusivity K, = a, Uj-"z!-™, which has 
been recently discussed by D. R. Davies (2). Also, since the diffusion 
equation can be solved for this particular problem using K, proportional 
to 2” and K, proportional to z!~”, this second method makes it possible 
to compare the effects of using these laws on the theoretical rates of 


evaporation. 


2. Diffusion of vapour from an infinite down-wind line source, and 
an empirical evaluation of the factor a, in the lateral diffusivity 
law K, a, Ui-"*z'— 

In order to solve the problem of diffusion from an array of infinite line 
sources, described in section 3, a solution is first needed for the problem of 
diffusion of matter from one infinite down-wind line source. A considera- 
tion of this latter problem also makes it possible to derive an empirical 
formula for a, in the lateral diffusivity power law by arriving at a solution 
In tWwO ways 

Professor O. G. Sutton has shown (4) that a continuous point source 
solution may be obtained by an application of his statistical theory to 
three dimensions, assuming wind speed constant with height. Despite this 
assumption, the solution gives the correct functional dependence on dis- 
tance from the source and also leads, in adiabatic conditions, to a close 
approximation to the absolute magnitude of height and width of the smoke 
cloud generated from a point source, if the kinematic viscosity of the air 
is replaced by the new concept of macro-viscosity. Sutton’s formula, using 
his notation, is 

2 f N 29,2 / (12 _|_ +2/(2)d 9 
x Og gan oe —2 2(y?/C2-+-2?/C2)}. (2.1) 
Suppose now a line source, of constant strength Q per unit length, 


extends from a point A on the x-axis down-wind a distance 2, to a point B. 
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Then the concentration of vapour at a point (y,z) in the plane normal to 
AB at B is given by 

2¢ q ; 

2@ | a"—ex p{ —ax"-2(y?/C2 +-22/C2)} da. (2.2) 


Y/) = a4 
C,%, © 


my 
0 


This may be written in the form 


22 [pfl=m)\_ pfy*/Cp+24/C2) 1—m)] - 
| lh 3 | 


x 7(2—n)C, Ca, 2 


= 6 (2.3 
(y? C3+-2? C2)a m)/2 \ ) 


x 
If the point A is allowed to recede to infinity in the up-wind direction, 
the solution for an infinite line source follows, 

— — 2QT{(l—m)/2} (2.4) 

/ Y 1/6 os : Y 242.(1—m)/2° inna 

aCmC,(2—n)i,{y?+ (C,/C,)222}0-m 

An alternative solution is now obtained by solving the appropriate 
diffusion equation, subject to suitable boundary conditions. Using the 


notation of previous work, the diffusion equation is written 


u ox — ; Kk, ox =e . K.2) (2.5) 
ex = by\ “Oy)  éz\ ez 


For an infinite line source, parallel to the mean wind direction, orien- 
tated along the x-axis, or for a semi-infinite line source, when z is large, 


the first solution shows that éy/éx = 0, and equation (2.5) becomes 


0 (Kv) i (K, | (2.6) 
oy cy O2z\ “co 


A particular solution of equation (2.6) is needed which fits the boundary 
conditions. The power laws found in recent work (2) to describe adequately 
the vertical and lateral diffusivities are written 

a, = 6,07" (2.7) 
and K, = a, Uj-*2'-*. (2.8) 
The application of these laws leads to good agreement with experiment for 
the decay of peak concentration variation from a continuous point source. 
Substitution of equations (2.7) and (2.8) into equation (2.6) leads to the 
equation 


0 = LE (z-mX)\ 4 | (= mX) (2.9) 
” Oy Oy} ez\ oz 
where a,/a, = f. (2.10) 


Writing y = »,/f, equation (2.9) becomes 


a2 ., a2 ¢ 
at WE x, m) Ox _ 0. (2.11) 


On* © O28 * 2 © 
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Now introduce polar coordinates given by 7 = rcos@ and z = rsin@, 
then equation (2.11) transforms to 
O*y Cx . oy CY 
p2o X (2 m)yr X14 X14 (1—m)cot aX = 0. (2.12) 
or or 70° o6 


A solution which is independent of @ is now sought in order to satisfy the 
boundary conditions of no mass transport across the plane z = 0, apart 
from the diffusion from the line source along y = 0,z = 0. If the equation 
is separated into normal solutions in r and @, it leads to Legendre functions 


which do not fit these boundary conditions. The equation in r is 


d*y l 
r=" X 1 (2 mr X = 0. (2.13) 
dr? dr 
The solution which tends to 0 as r tends to infinity is 
x = Ajri-m, (2.14) 


where A is a constant determined by the conditions at the source. Then, 
in the original y and z coordinates, the form for x is 


x Af 1—m)/2 (y” } 2) -m)/2 (2.15) 


The constant A may be determined by considering the vertical flux of 


vapour across some plane 2 m1, Say, at a lateral distance y, i.e. 


(— K, éx/0z),.,,- 


Then integration with regard to y from negative to positive infinity will 
give the total flux across a strip, parallel to the ground, unit length down- 
wind and infinite across wind. This is easily seen to be independent of z, 
and must be equated to Y, the strength of the source. The resulting 
equation leads to an expression for A in terms of Q, m, f, a,, and U,, the 
mean wind speed at a given reference height. The result needed here is 
that y is proportional to 


+2\(1—m)/2 
y/Gz)2"5 : 


L/jy*+(a 

It is reasonable to suppose that the two methods lead to the same 
approximate functional dependence on y and z, and this is seen to be the 
case, ifa,,/a, from equation (2.15) is identical with C?/C? from equation (2.4). 
From the definitions of C!, and C, given by Sutton (4) this identity leads 
to the equation 
(2.16) 
The ratio v’®/w’? is known to have a mean value of about 1-8 over the 


lowest 5 m. of the atmosphere (4), in adiabatic conditions, and this leads 
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to a value for a, of about 1-6a,. The factor a, may be evaluated directly 
from meteorological measurement and a, then derived. It has not yet 
been possible to verify this empirical formulation of a, for evaporation 
problems using 


4 Tl—-n+l—m 
K, a, l ‘aye, 


but it is an experimental fact that, in adiabatic conditions over downland, 
cloud height and semi-width at a distance of 100 m. down-wind of a con- 
tinuous point source are 10 m. and 17-5 m. respectively; this indicates that 
the lateral and vertical diffusivities are approximately in the ratio of 1-75: 1 
over the lowest 10 m. of the atmosphere. Hence, in adiabatic conditions, 
the application of equation (2.16) leading to a value for a,,/a, of 1-6 appears 
appears in all 


to be in fair agreement with observation. This factor, a,, 


evaporation calculations involving the effects of lateral diffusion. 


3. The problem of turbulent diffusion of vapour from a plane, 
saturated, evaporating area of constant, finite width, when @y/éx 
is assumed to be zero. Solution by the method of ‘sources’ 

It has been assumed in evaporation theory that, for wind speeds greater 
than a certain low critical value, the integrated rate of evaporation is 
simply proportional to width to the power unity. However, visual observa- 
tions of a saturated, rectangular-shaped piece of filter-paper placed tan- 
gentially in the air stream of a wind tunnel show very clearly that the local 
rate of evaporation near the leading and longitudinal edges is considerably 
higher than it is elsewhere. The theoretical results obtained in the case 
of a parabolic area (1) also point to this effect. It thus appears likely that, 
as the edge effect will be sharper for areas of small widths than it will be 
for larger widths, the integrated rate of evaporation per unit length of the 
evaporating strip will be proportional to the width of the strip to a power 
less than unity, for any wind speed. 

An array of elemental infinite line sources is now considered, arranged 
in the plane of the ground so that each is parallel to the x-axis and con- 
tained in the region 0 < y < w,z = 0. The strengths, Q, of all line sources 
are taken to be constant along the length; the concentration due to the 
particular line source lying, say, along the a-axis, is given by equation (2.4). 
Suppose now Q is a function of y (distance across wind) and choose Q(y) 
so that y is constant in the region z = 0,0 < y < w. The integral equation 
obtained for Q(y) is closely allied to the Abel type; a solution can be ob- 
tained numerically. 

Consider a point P (see Fig. 1) on the area (z = 0), composed of the 
system of line sources, and at a distance 6, say, from the edge y = 0. 
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Now the concentration at P due to a line source at a distance y from the 


z-axis is obtained by putting z = 0 in equation (2.4). The value is 
Xx AQ(y)/(@—y)™. (3.1) 
,{1—m etd " 
where A 2] | = ) |x yp C(2—n)it. 
ul 
Or Wo anna nn nnn nnn ey 
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Consequently the total contribution to y at P due to the whole array of 
line sources from y 0 to y wis 


H uw—Oé 
| P (y) dy ; UY) dy 
d A y)} m | (6 y)} m\i* 
0 0 
This leads us to an integral equation for Q(y), 
w—@é 
V y dy [ VU(y) dy Xo 
(é . 


(A y)} m } y)} m A 


(3.2) 


Equation (3.2) is closely allied to Abel’s equation, but an analytic 
solution does not appear to be readily obtainable. However, a solution 
for Y(y) can be deduced from the results of the following section, and 
tested here numerically in particular cases. 


Assume O(y) B/y™?(w—y)"?, (3.3) 
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where B is a constant to be determined by substitution into the integral 
equation. In a particular case given by w = 20 units, m = 1/7, it has been 
checked numerically to three significant figures that the function (3.3) leads 
to the same value of y at y = 0, 3, 5, 7, and 10, as expected; the value 
of x,)/AB was 14-1 in each case. The concentration x at any point above 
the system of line sources can be expressed formally as an integral once 
the value of the constant B has been evaluated numerically from the 
integral equation in each particular case. The interest in this paper lies 
in the dependence of integrated rate of evaporation EF, per unit length, on 
the width of the area. In equation (3.3) write y’ = y—w/2, then Q is seen 
to be proportional to 1/{(w/2)?—y’*}"?, and E is proportional to 

La) 

| {(de)?—y’2}-m? dy’. 

w/2 
By writing y’ = (wsin@)/2, E£ is easily seen to be proportional to w1-™. 

Since we have chosen Q(y) in such a way that x is constant over the area, 

the conditions are similar to those over a saturated, plane, evaporating 
surface of finite width, bounded by straight edges parallel to the mean wind 
direction, and it is in these conditions the result is obtained that EZ is 
proportional to w1-™, 


4. Second method of obtaining a solution to the problem of 
turbulent diffusion of vapour from a plane, saturated, evaporating 
area of constant finite width when ¢y/éx is assumed to be zero 
In the problem outlined in section 3 the conditions over the area com- 

posed of line sources are (@) y = y, on the area, (b) éy/ex = 0 everywhere, 

and (c) no mass transport across the plane z = 0, outside the area. Use 
is now made of the power laws for K,, and K,, given by equations (2.7) and 

(2.8). On substitution of (2.7) and (2.8) into equation (2.6), the appropriate 

transport equation is written, as in section 2, in the form 


0, (4.1) 


where 7 = (f)-4y and f = a,,/a,. 
In order to obtain suitable solutions of this equation, elliptic coordinates 
are introduced, given by 
n+iz = ccosh(B+76), (4.2) 
where 2c is the total width of the strip in (,z) coordinates. Equation 
(4.1) transforms to 


apt opt m)eoth B25 + (1—m)cot =F = 10. (4.3) 
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Appropriate boundary conditions are 
(i) x/X9 > 1 as B> 0, x, being the saturated value for x at the sur- 
face of the lake, 
(ii) y>0 as B>o 
(ili) 0x/e8 > 0 at O= 0, z. 
A solution of equation (4.3) exists, independent of @ (and consequently 


fulfilling condition (iii)), satisfying the equation 


qx 1 (1—m)coth B2X = 0. (4.4) 
dp dg 

It is of interest to note here that normal solutions of equation (4.3) in 
6 and 8 coordinates are the Legendre P and Q functions, but these func- 
tions do not fit the boundary conditions as they are usually formulated 
for evaporation problems. Equation (4.4) is easily integrated to fit the 
boundary conditions, leading to the equation for x, 


| (sinh 6)"—1 dd 
} 0 ; 4.5 
ane I(m) aaa 
where I(m) ( (sinh d)"-! dd. 


0 
These integrals are convergent in the neighbourhood of B = 0 and at 
5 infinity for m Q. 

[t is of interest to note that if we put » = 0 in the formula previously 
obtained by D. R. Davies (1, equation 3.1) for the distribution of vapour 
over a parabolic area, the solution for the problem immediately reduces 
to the form of equation (4.5) of the present paper. 

Equation (4.5) may be thrown into a form more suitable for computation 
as follows 


Write sinh dé = tan uw, (4.6) 


and on substituting, the integral J(m) becomes 


I(m) | (cos u)~™(sin u)™—-! du. (4.7) 
0 
This may be expressed in terms of [ functions, by using a well-known 
result (5, p. 256), in the form 
I(m) r{4(1- m)}T(4m)/21 (4). 
The numerator in equation (4.5) may be computed numerically to a 
sufficient degree of accuracy by using Simpson’s rule: this was the method 
used satisfactorily in the first instance for a similar integral occurring in 
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the solution of the parabolic area problem (1). The mathematical form 
of the solution found in the present problem is very similar to that found 
in the previous work, but it has not been found possible to transform it 
into a convenient series of well-known functions. 

A very simple expression may now be derived to give the rate of evapora- 
tion per unit length of the lake in terms of [ functions. From equations 
(4.2) and (4.5) it is easily calculated that, as z > 0, 


K, 6x /0z > x,a, U{-"/e™ I(m)(sin 0)". (4.8) 
Now at 8 = 0, 7 = ccos@ and hence 
\dy| = (f )tesin é dé. 


Consequently the net rate of evaporation per unit length becomes 
hor 

) Tl-npl—mfs 

. 2y,a, cl-mf? : 

E Xone 1 (sin 6)!-™ dé 

I(m) 


“ 


Xo 4, UZ-"cl-™T'(3) (1 —4m) ft (4.9) 
I(m)P(3—4m) “4 


In this expression 2c refers to the width of the strip in coordinates. In 
the original y-coordinates the formula becomes 
2ry,a, U1-"f @T(1—}m)w!-™ 


= (4.10) 


E — 
T(3m)0(8—4m)0(4—4m) 


where w represents in appropriate units the half-width of the lake, i.e. Z 
is proportional to wi-™. 

[t is interesting to note in this section that, if we make use of the power 
law for lateral eddy diffusivity, K, = a, Uj-"z™ (discussed previously by 
D. R. Davies (2)), instead of that given by equation (2.8), the above 
analysis applies after an initial transformation of coordinates; the details 
are described in the Appendix to this paper. One result of this calculation 
is that the theoretical dependence of rate of evaporation on width is given 
by # proportional to w/@+2™; in adiabatic conditions this becomes w8 
compared with w®**, using the z!-” law. In the Appendix the theoretical 
formulae for rates of evaporation are worked out numerically and com- 
pared for the two laws, assuming a particular set of conditions, and they 
are found to agree very closely. 


5. Experimental arrangements and results 

Some experiments have been carried out by the authors in an open- 
circuit wind tunnel to test whether the theoretical law for the dependence 
of evaporation on width obtained in sections 3 and 4 is approximately true 
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for plane, rectangu!ar, saturated surfaces of finite length in wind-tunnel 
conditions. 

Visual observations of a saturated, rectangular piece of filter-paper, 
stretched on a flat glass plate which was held by a faired holder near the 
floor of the working section of the tunnel, showed very clearly that, as 
the area begins to dry, the rate of evaporation was very much greater in the 
region near the leading and longitudinal edges than it was elsewhere. This 
is as we should expect for two reasons: (a) the vapour is transported from 
the edges in directions varying from those parallel to the glass plate to 
those normal to the plate; (b) the concentration is smaller over the edges 
than it is in the centre, and, incidentally, is smaller than it would be if 
the area was infinite across wind, and consequently there is less ‘resistance 
to evaporation’ in this region than in the centre. 

Quantitative experiments have also been made to test the dependence 
of total evaporation from an area on the width of the area. Exhaustive 
tests had previously been made by R. W. Powell and E. Griffiths (6) to 
find the dependence of evaporation on length and width. They found a 
law of the type given in sections 3 and 4 of this paper for dependence 
on width in low wind speeds, but the index approached unity as they 
increased their wind speed. For wind speeds greater than about 300 cm./sec. 
they conclude that evaporation is exactly proportional to width and so 
lateral diffusion may be neglected. However, the evaporating area was 
mounted, during the experiments by Powell and Griffiths, in the tunnel 
midway between the floor and ceiling in a special holder fitted with a 
streamlined nose. Experiments by Pasquill (7) on this type of suspension 
indicated strongly that over the leading portion of the evaporation plate 
the boundary layer is laminar, changing gradually through a transition 
zone to a turbulent state. The curves of velocity and vapour-pressure 
dependence on height given by Powell and Griffiths (6, Fig. 10) indicate 
strongly that a transition from laminar to turbulent flow does indeed take 
place towards the leeward edge of their plate. In the theory described in 
the present paper the presence of a fully developed turbulent boundary 
layer is assumed, and to test the theory this condition must obtain over 
the evaporation plate. 

The experimental technique followed is that used by Pasquill in experi- 
ments described fully in his paper (7). Rates of evaporation, air tempera- 
ture, and wind speed in the free stream of the working section were 
measured as in his experiments. Recorded variations in the air tempera- 
ture over a run of 3 or 4 hours were very small, and corrections were made 
for any deviations from a chosen standard, following the method outlined 
by Pasquill. Rectangular shaped pieces of good filter-paper (Whatman 
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No. 1) were saturated with aniline and stretched on a plane sheet of glass 
(8 in. square) which was placed centrally, on the floor of the working 
section of the tunnel (2 ft. 6 in. section), in a faired holder positioned 
7 ft. from the inlet end of the tunnel. In this way a fully developed 
turbulent boundary layer was formed on the plate. Measurements of rate 
of evaporation, #, in grammes per minute, were made from rectangular 
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Fic, 2(i).‘ Variation of log E with log w, when 

length of area 5-0 in., mean air speed 

1,330 ft./min., air temperature = 16-0° C., 

giving FE oc w84, EF denotes rate of evapora- 

tion in g./min., and w denotes the width of 

the area in in, 
areas by the ‘gravimetric’ method used by Pasquill. The indicating balance 
used was shown by direct test to be reliable to 0-005 g. From initial tests 
it was possible to calculate the time, in a given set of circumstances, during 
which the rate of loss remained constant before the ‘drying period’ set in. 
Losses of about 1-8 g. were recorded from the larger areas used, and for 
the smallest areas losses were about 0-30 g. The mean air speed in the 
tunnel was kept constant during each run and its value recorded by means 
of a calibrated air meter positioned in the centre of the tunnel, downwind 
of the evaporating area. The mean air speeds used were greater than those 
considered to be critical by Powell and Griffiths. The wind velocity 
‘profile’ near the evaporation plate was investigated by means of a small 
standard pitot static tube, these observations yielding an m value of 0-15 
An appropriate value for (1—m) is thus 0-85. 

For a constant length of area and a constant mean air speed log E was 
plotted against log w, w being the width of the area in inches. The results 
obtained are given in Fig. 2 in the form of four sets of points, each set 
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corresponding to a given length and a given mean air speed; the various 
conditions tested are described in the figure. A straight line was drawn 
through each set of points, and in the case given in Fig. 2 (i), where twenty- 
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Fic. 2 (ii). Variation of log E with log w, when length 

of area 6-75 in., mean air speed 1,330 ft./min., 

air temperature 19-3° C., giving FE o w®®*, E is the 

rate of evaporation for the area in g./min. and w is the 
width in in. 


Log w 








Fic. 2(iii). Variation of log FE with log w, 


when length of area 3-5 in., mean 
air speed 1,330 ft./min., air tempera- 
ture 16-5° C., giving Eo w*®?, 


five points were plotted, a line of regression was computed. The slopes of 
the straight lines are found to be 0-84, 0:80, 0-80, and 0-81 in the cireum- 
stances given in Figs. 2 (i), 2 (ii), 2 (iii), and 2 (iv), respectively. 


The theoretical law of dependence of evaporation on width (i.e. 
sl m 
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w®'®5 in the circumstances tested) is thus in reasonable agreement 
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with experimental values (between 0-80 and 0-84) in these wind-tunnel 
conditions. 
40g w 
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of 


area 6-75 in., mean air speed 720 ft./min., air tempera- 


Fic. 2(iv). Variation of log # with logw, when length 


ture 17:0° C., giving EF oc w®!, 


6. Conclusions 

The mathematical form of the final theoretical expression for rate of 
evaporation per unit length of an area gives the kind of dependence on the 
various entities to be expected. The theoretical result that H is propor- 
tional to w'-™ is of particular interest, and experimental results recorded 
in wind-tunnel conditions up to 700 em./sec. indicate that this kind of 
relationship is also approximately true for rectangular areas, even when 
the length of the area is not large compared with the width. It is interest- 
ing in this connexion to note that an examination of the solution obtained 
for the parabolic area problem shows a similar result in that the rate of 
evaporation for an elemental strip of the parabola, in the cross-wind 
direction, is proportional to width of the strip to the power j, if we put 
m = }, in adiabatic conditions. 

[It is not claimed in this paper that the absolute values for rate of 
evaporation given by equation (4.10) is even approximately true in the 
region of a rectangular area where the effect of the leading edge is important: 
thus a calculation of equation (4.10), assuming the conditions given in 
Fig. 2 (i), leads to a result which is only 62 per cent. of the observed value. 
The importance of taking an index (1—m) instead of 1 in the width law 
lies in its application to the problem of using observed values of rate of 
evaporation from small areas to predict rates of evaporation from large 
areas. This problem is complicated by many factors, but the law of 
dependence on width obtained in this paper is certainly one which must 








ar WN = 


i 


al 


the 
0r- 
ded 
| of 
hen 
est 


ned 


ind 


put 


nust 














EFFECT OF FINITE WIDTH OF AREA ON RATE OF EVAPORATION 479 


be allowed for, as the effect of using an index of 0-84 say, instead of 1-0, 
is evidently very considerable for large widths. 

Finally, it is noted that the method of ‘sources’ followed in section 2 
may have application in other problems where edge effects are important. 
In the case of evaporation from a parabolic area, for example, the distribu- 
tion of rate of evaporation, and consequently of source strength is already 
known theoretically (1) and the method of ‘sources’ then leads us imme- 
diately to a formal solution, expressed in terms of double integrals, for the 
distribution of vapour down-wind of the area. 
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APPENDIX 
A COMPARISON OF THE EFFECTS OF USING THE 2” AND 2!-™ LATERAL 
DirFusIvITy LAWS ON THE THEORETICAL DIFFUSION OF VAPOUR FROM 
A SATURATED AREA OF FINITE WIDTH, WHEN @x/0x IS ASSUMED TO BE 


ZERO 
Assuming that the lateral eddy diffusivity may be written 
K a, U8", (1) 
and writing a,,/a f’ and y jf’, the appropriate transport equation becomes 
(2m &X) . S (-1-mEX) _ 9, (2) 
CY cy Cz C2 


Now writ =m+} (m-+4)e, (3) 
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and equation (2) becomes 
ay ax  1/(2m+1) ex 0 4 
oy” Toyz © v ov : (4) 
This equation is identical in form with the basic equation (4.1) arrived at in 
solving this problem using the z!—” law, if the expression (1—m) in the z1—™ case 
is replaced by 1/(2m-+-1). By making the transformation into elliptic coordinates 
y’ +iv ccosh(B+76), (5) 
the ensuing relation between B and y/yp is obtained, following the method previously 
used in section 4, of 8 
| (sinhd)-Vem+) dg 
0 
sn (6) 
Xo I’(m) , 
where I’(m) D{1/2(2m-+ 1)}P{m/(2m+ 1)}/2T (4). (7) 
If the details of the solution based on the z!~” law are pursued, the following final 
result is arrived at, using the same notation, for the net rate of evaporation per unit 
length of the area, 
2X0 a, uw: os A at (2M+1)q4y1 (2m+1) (9m | })} (2m+1) 


- Tim (2m + 1)}1{(4m+3)/(4m-+ 2)} 


: (8) 
This result may be obtained from equation (4.10) on replacing (1—m) by 1/(2m+1) 
and multiplying by a factor of (m+ 4)/@™+)_ There is little difference in the depen- 
dence of E on width of area in the two cases: in adiabatic conditions E is proportional 
to w®®§ and to w®’® for the z!~” and 2” laws respectively. Denoting the solutions 
using the z!-™ and z™ laws by Z,_,, and E,, respectively, the ratio of the two solutions 
is found to be ‘ 
( f’ym 2m-+ 1) (qyy)1/(2m-+ NT km (3 — km)(m te. ) (2m+1) 


i ci {1 m \_(4m+3 , : (9) 
a Cf )Pm(w)-mr | =)r( ) 





2m+1)" \4m+2 

This ratio has been evaluated, assuming that adiabatic conditions are represented 
by m + and taking w = 500cm. The factor f, following the indications of section 2 
of the present paper, is taken to be 1-6. The factor f’ is found to have the value 44-0, 
assuming that a, and a, are given sufficiently well in the 2” case by equations (2.3) 
and (2.5) of previous work by D. R. Davies (1)—they certainly lead to good agree- 
ment between theory and experiment at the downwind edge of a 10-yard long 
parabolic evaporating strip in the conditions tested. The resulting numerical value 
for E,/Ey_m is 1-005. 


The z!-™ law of lateral diffusivity leads to the correct continuous point source 


m 


solution for decay of peak concentration, and therefore it is reasonable to assume 
that it will be satisfactory for calculations of rates of evaporation. Hence the result 
that the z” law gives a numerical value for rate of evaporation, in these particular 
circumstances, within } per cent. of that obtained using the z!-™ law, provides 
further support for the use of the 2” lateral diffusivity law in evaporation problems. 
In wind-tunnel conditions it is known that this law leads to theoretical results in 
reasonable agreement with observations for transfer of vapour from small, parabolic, 
saturated areas; and in the atmosphere, if a, is chosen so that K, = K, at 200 cm. 
above the surface, the 2” law has been successful in predicting the distribution of 
vapour above parabolic areas. The advantage of using this z” law to represent 
lateral diffusivity in evaporation problems lies in the fact that the resulting diffusion 
equation is much more amenable to mathematical treatment than the equation 
following the use of the z1-” law. 
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